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FORE WORD 

AN Industrial Computation Seminar, sponsored 
jL JL by the International Business Machines Corpora- 
tion, was held in the IBM Department of Education, 
Endicott, New York, from September 25 to September 
29, 1950. The ninety research engineers and scientists 
who participated in this Seminar met to discuss the 
fundamental computational methods which are appli- 
cable in a wide variety of research problems. Particular 
attention was drawn to computational techniques re- 
cently developed in the fields of chemistry and petroleum. 
The International Business Machines Corporation wishes 
to express its appreciation to all who participated in 
this Seminar. 
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AS I LOOK over the list of occupations of the members 
of this Seminar, I am impressed by the wide range of fields 
represented — engineering, physics, chemistry, accounting, 
even astronomy. It is interesting to note how experience in 
one field can influence seemingly unrelated activities in 
other fields. I might illustrate by a trivial example of a 
procedure mentioned here in connection with an accounting 
problem which also occurred recently in an astronomical 
problem at the Watson Laboratory. It was mentioned that 
it is frequently more convenient to produce the calendar 
date on the accounting machine than to copy it from the 
calendar. In our problem we required calendar dates at 
forty-day intervals from 1653 to 2060, taking into account 
the complicated leap-year rules of our calendar; the list 
was prepared by a single run on the IBM Type 602- A Cal- 
culating Punch. 

The close relationship between apparently unrelated 
things is not new in science; it is necessary to look at the 
picture in the proper perspective to see the relationship. 
Let us consider the ancient astronomer who was intrigued 
by the small spots of light in the sky, called planets. He 
spent many hours measuring and recording their positions, 
and I am sure that his contemporaries could not see how 
these activities would ever put food into the mouths of men. 
His contemporaries could not see enough of the picture. 
We now know how the study of such planetary observa- 
tions has led not only to our understanding of the motions 
of the planets, but also to our knowledge of the funda- 
mental principles of mechanics, the basis of all mechanical 
design. From this vantage point, when we see a farmer 
riding a properly designed tractor and pulling carefully 
designed implements, we realize that the early astronomer 
has done more to feed the multitudes than all of his con- 
temporaries. 

The computing profession has always incorporated math- 
ematical and mechanical techniques; benefits from one field 
of science have been carried into other fields. At the time 
of the early astronomical observations computing was being 
done on a considerable scale, and from that time forward 
man has tried to develop computing aids. The mathema- 
tician and the scientist have tried to devise both mechanical 
aids and mathematical aids. The first astronomer needed 



trigonometric tables in order to make his computations. 
Later, the development of the logarithmic table greatly 
facilitated his arithmetic operations. The adding machine 
was invented by Pascal in 1642 and the desk calculator by 
Leibnitz in 1693. Thus, we have the invention by scientists 
of these two tools that were greatly needed by scientists; 
yet they were of little use to science for over two centuries. 
Although the desk calculator of Leibnitz was, in principle, 
our present-day machine, two centuries were required to 
'develop it into a generally useful implement. 

There are two reasons for this long delay: one is mechan- 
ical, and the other is mental. It is a long, hard pull from 
the gleam in the inventor's eye, or even from the first 
model, to the point where a scientist can use a device as a 
help and a tool, rather than as a problem in itself. It is not 
difficult for the inventor to make his model work under his 
own benevolent criticism, but to have it developed to the 
point where it is accurate, fool-proof, and efficient is an- 
other matter. Of course, things did not proceed as rapidly 
in those days as they do now, but one should not overlook 
the great number of technical developments necessary to 
fill in the details that make a complicated machine work. 

On the mental side there was the fact that people had 
been trained in the use of logarithms, and computational 
work had been organized for the use of logarithms. If you 
look through many books in applied mathematics of that 
era, you will find that large portions of many of them were 
devoted to the conversion of basic formulae into a form 
suitable for the effective use of logarithms. Then, too, there 
was the matter of tables. To replace logarithmic tables with 
natural tables required some time. This seems like a mod- 
ern age, yet I am not an octogenarian and I can remember 
the dying gasp of the logarithmic table as the standard 
method of computation. I have seen the desk calculator be- 
come a necessary instrument for every scientist who is 
doing quantitative work. They are now so efficient and so 
reliable that the scientist has merely to insert the proper 
number, push the proper control key to perform his desired 
operation, and read the results. 

In medicine there is a well-known phenomenon. While 
smallpox and diphtheria were decimating mankind, little 
attention was paid to some of the lesser diseases, but now 
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that smallpox and diphtheria have been brought under con- 
trol the lesser diseases are considered very important. Simi- 
larly in computation, when multiplication and division were 
performed longhand, or with the aid of logarithms, the 
computer did not worry about the associated clerical oper- 
ations. The advent of the desk calculator, however, enabled 
him to undertake larger problems, and he has become pain- 
fully aware of the details of reading and writing data and 
of initiating the proper control operations. The need for 
automatic handling of data and instructions becomes impor- 
tant for further progress. 

Here again, we have a situation similar to that of the 
development of the desk calculator. In 1893, as you have 
heard, the punched card was introduced as a means of read- 
ing data and instructions into a machine, but many years 
elapsed before the method became accepted as a regular 
part of scientific research. This delay also has been due to 
the necessity of both mechanical and mental development. 
When Mr. Watson became president of the International 
Business Machines Corporation in 1914, he immediately 
organized a development program to make the machines 
more versatile. From 1914 to 1930 there was extensive de- 
velopment of the various functions of the machines such as 
reading cards, sorting, printing, adding, subtracting, multi- 
plying and recording. During the following two decades 
further remarkable technical development in the machines 
occurred, but in my opinion the mental development among 
the scientists and engineers has been more striking. During 
this time there has come the general realization that the 
punched card had already provided the means of automati- 
cally handling scientific and engineering data. The recent 
introduction of the electronic circuit, which has greatly in- 
creased the speed of some operations, has dramatized the 
automatic process, but the hundreds of successful automatic 
computing installations now in operation have their roots in 
the mental revolution of the past two decades. 

It is interesting to note that during this period Mr. 
Watson saw the importance of automatic computation to 
science more clearly than the scientist or the engineer, and 
the present widespread use of such facilities is due in large 
measure to his early efforts. These efforts included not only 
the rapid development of standard machines to make them 
more generally useful, but the development of many special 
devices for academic purposes. In 1928 he established the 
Columbia University Statistical Bureau for educational re- 
search, and soon after installed there a special statistical 
calculator. The operation of this machine was very striking 
even by today's standards. It would read data and limited 
operating instructions from the cards at the rate of nearly a 
million digits an hour; it would add 100 digits simultane- 
ously from the cards or from other parts of the machine 
according to a complicated program, and print the results 
at the rate of nearly one-half a million digits an hour. In 
1933 a second laboratory was established at Columbia capa- 



ble of handling general scientific calculations such as the 
solution of differential equations, and the reduction of ob- 
servational data. This laboratory was in full-time operation 
on basic research until the advent of the war when it was 
converted to military research. By 1940 a number of labo- 
ratories about the country were using standard punched 
card machines for technical computation. During the war 
new laboratories were quickly established in all phases of 
the defense effort, including atomic energy, aircraft design 
and construction, air and sea navigation, and many others 
that are still classified. In 1944 the IBM Automatic Se- 
quence Controlled Calculator was completed in Endicott 
and presented to Harvard University; this machine, known 
as Mark I, has since been in continuous service. In the 
same year two relay calculators were installed at the Aber- 
deen Proving Ground; these machines were more limited in 
capacity and flexibility but were about twenty times as fast 
as the sequence calculator. They are still the fastest relay 
calculators in operation. 

In 1945, a special table printing device was installed at 
the Naval Observatory which enabled the scientist to print 
mathematical tables from punched cards in a form suitable 
for direct reproduction by the printer. In the same year, the 
recording equipment for the great wind tunnel at California 
Institute of Technology was installed so that the observa- 
tional data could be recorded directly in cards without 
hand transcription. Then came the 603, the first commercial 
electronic calculator, which has been replaced by the more 
versatile IBM Type 604 Electronic Calculating Punch. The 
Selective Sequence Electronic Calculator, which was dedi- 
cated in January, 1948, provided electronic speed of opera- 
tion together with an internal storage capacity of half a 
million digits and completely automatic programming. The 
Card Programmed Electronic Calculator is the most recent 
addition. 

From such a wide variety of available punched card ma- 
chines- — sorters, accounting machines, collators, reproduc- 
ers, the 602, the 602-A, the 604, the CPC, and the SSEC— 
many of you are probably wondering which equipment you 
should use and why. Since much of the later part of the 
program deals with the more recent machines such as the 
604 and the CPC and their detailed application, I shall con- 
fine my attention to some of the more basic uses of the 
punched card and of the simpler punched card machines. 
I shall take my examples from some of the earlier work; 
it is interesting to note that in many applications these early 
techniques are still the most efficient in spite of all the ad- 
vances that have been made in design. Moreover, many of 
these early techniques clearly illustrate basic principles that 
can be readily applied in general. 

A question of basic importance in the application of the 
punched card is whether the calculation should be done 
sequentially or in parallel. In a computation with a desk 
calculator the problem also arises, and we shall illustrate it 
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with a simple calculation. Let us assume that we wish to 
evaluate the formula 

f{x) = a -{- bx -\- ex- + d sin(a + hx + cx"^) . 
Ordinarily we should write the formula across the top of 
the sheet and assign successive columns for the intermediate 
results. Successive rows would be assigned for successive 
values of x. We could perform the computation by evaluat- 
ing the formula completely for the first value of x, including 
the use of the sine table. In this case the work would be 
completed a row at a time, or sequentially. On the other 
hand, we could perform the first operation, say b • x, for 
all values of x before proceeding to the next operation. This 
procedure by columns is termed parallel computation. The 
experienced hand computer knows that he can perform 
operations such as table lookup most efficiently by the 
parallel method. The user of the desk calculator, however, 
can frequently avoid the recording and reading of inter- 
nnediate results by using the sequential method of operation, 
for example, the formation oi a -{- bx -\- cx'^. Thus, he might 
combine the sequential and the parallel methods. In the use 
of the punched cards, parallel operations become more im- 
portant because of the ease of reading and recording inter- 
mediate results in the cards. 

The ability to store numbers temporarily in the machine 
largely determines the length of the sequence that can be 
handled before the recording of intermediate results. The 
602-A, for example, can store a dozen numbers while the 
SSEC can handle 20,000. Perhaps we should all like to use 
a, machine that can handle a million digits, but economic 
considerations must be taken into account. It must be re- 
nnembered that a small box of cards costing about ten dol- 
lars has more storage capacity than a million-dollar calcu- 
lator. The desirable amount of storage in the machine is 
thus an economic matter, and the amount of storage that is 
rieeded depends to a considerable extent on the nature of 
the problem. Where parallel computation can be employed, 
a.nd a very large portion of technical work can be so han- 
dled, the small machine with card storage is indicated. 

The length of sequence that a given machine will handle 
depends on the facilities for handling instructions as well 
ais on the storage capacity. On the 602-A it is easy to wire, 
say, ten, twenty or thirty successive operations; many more 
than that are possible, but the wiring becomes somewhat 
troublesome. The CPC has greater storage capacity than 
the 602-A and also the ability to read operating instructions 
from punched cards. It can therefore handle sequences of 
any length, limited only by the storage capacity. The SSEC 
has a tremendous storage capacity and facilities for the 
most intricate programs. Not only can instructions be read 
from cards, but the entire storage capacity of the machine 
can be used for manipulating the operating instructions. 

Efificient use of the punched card also requires that full 
advantage be taken of the facilities for rearranging data. 



The simple sorter can rearrange cards at the rate of 500 or 
600 a minute. When you consider that each card carries 
eighty digits, it is easy to see that in a day you can rearrange 
and reassociate millions of digits. Previous to the punched 
card there was no facility for automatically rearranging 
data; the computer wished to avoid such manipulation 
rather than encourage it. In adopting a new medium for 
computation it is desirable to profit from past experience, 
but not to limit one's thinking by the limitations of the past. 

To illustrate new possibilities introduced by the ability to 
rearrange data, I shall describe a large-scale operation car- 
ried out in 1928. The results of a survey contained 50,000 
cases with a dozen variables. These data were punched in 
50,000 cards, and in a week or two all correlations between 
all of the variables were computed. Those of you who are 
familiar with correlation analysis know that this work in- 
volved the formation of the sums of millions of products. 
There was no multiplying machine available, but by sorting 
the cards and adding from them in the proper order, we 
obtained the required results and printed them automati- 
cally. A small mental calculation shows the efficiency of 
this method. By feeding eighty-digit cards at the rate of 
9,000 cards an hour it is possible to accumulate partial 
products at the rate of 720,000 digits an hour. It is neces- 
sary to make a run for each digit of the multiplier and to 
allow space in the counters for accumulation. Allowing for 
three-digit factors it is easy to obtain 40,000 products an 
hour, a considerable achievement even in modern terms. 
This method, known as progressive digiting, or similar 
techniques continue to recur as an effective method in 
many problems including harmonic analysis and the forma- 
tion of normal equations. Since the method requires only a 
key punch, a sorter, and an accounting machine, it has been 
applied where only simple accounting installations are 
available. 

An important change in emphasis brought about by the 
ability to rearrange data is in the use of mathematical tables. 
In hand computation the use of tables is laborious and is 
frequently avoided. The punched card, however, makes the 
consulting of tables one of the most efficient operations. The 
sorter will rearrange cards, and each card may contain 
eighty columns of tabular data. For example, a single card 
can carry an angle, the sine, the cosine, the tangent and any 
other data associated with the angle. The sorter will put the 
cards in the proper place, and the reproducer will transfer 
the data into the computing cards at the rate of 6,000 cards 
an hour. This means that with this simple equipment you 
can look up half a million digits of tabular data in an hour. 
This ability to find millions of digits of tabular data in a 
day with the simplest equipment and a few cards is so fun- 
damental that the scientist must re-examine his whole con- 
cept of computation. He no longer says "How can I avoid 
table lookup ?" but "How can I use more tables ?" 
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Before leaving the subject of tables, I might mention one 
example of their use. In 1940 the U. S. Naval Observatory 
in Washington had the task of producing the American Air 
Almanac for the use of navigators at sea and in the air. It 
contains 730 pages a year, and each page contains several 
thousand digits. Since the lives of the crew depend upon it, 
the highest accuracy is demanded. The entire computation 
was performed automatically on standard punched card 
machines. Three months elapsed from the time the book 
was decided upon until the machines were delivered and in 
operation. The necessary planning and preparation of basic 
data had been done meanwhile. The first volume was pre- 
pared and checked, published by the Government Printing 
Office and distributed to the fleet before the first of the fol- 
lowing year. For the benefit of those who instinctively think 
of the most recent and most powerful equipment for a given 
problem, the work here described was performed with only 
the key punch, the sorter, the accounting machine, and the 
reproducer-summary punch. The multiplier was available, 
but it was not used because the other method was more 
efficient. The whole volume was produced by continuous 
summation in parallel columns involving the sexagesimal 
system. The summations were recorded for a ten-minute 
interval, and the data to be summed were obtained from 
special punched card tables prepared for the purpose. The 
standard accounting machine performed the accumulation 
in the sexagesimal system by means of fractional counters 
that carry on 6 instead of 10 (these fractional counters are 
regularly available) . After the initial volumes had been pre- 
pared, a special table printing device was installed to tran- 
scribe the resulting data from the cards in a form suitable 
for reproduction by means of line cuts and electroplates. 
In the twelve years that the publication has been produced, 
there have been computed and published about twenty-five 
million digits, which have been examined and used by thou- 
sands of navigators, and as yet not a single error has been 
reported. 

Another example of the use of simple punched card 
equipment dates back to 1935; this is a large-scale reduc- 
tion of observational data, which should interest many here 
who are associated with organizations that regularly handle 
large quantities of industrial, engineering, or scientific data. 
Here, again, we have an example of a procedure that applies 
equally well in any field; the punched card machines cannot 
tell the difference between data concerning a chemical com- 
pound or the wing of an airplane. The program handled in 
1935 concerned the measurement of astronomical photo- 
graphs. Each plate contained several hundred black spots 
or star images. Since the blackness and size of the image 
depend upon the brightness of the star, the brightness may 
be determined by the measurement of the amount of light 
obscured by the image in a photometer. The total program 
involved several hundred plates with several hundred thou- 
sand images to be measured- The problems involved were 



those common to all measurements: identifying the object 
to be measured, recording the instrumental reading, cali- 
brating the instrument, applying instrumental corrections, 
performing mathematical transformations, combining re- 
sults, discussing the errors statistically, and recording the 
results for publication. The punched card method permits 
all these operations to be carried out more thoroughly than 
before, without error, and with a minimum of drudgery. 
The identification of the image at the instrument was made 
with the aid of scale settings computed for each image on 
the card; these settings were printed on the card by the 
interpreter. The observer placed the card in a punch and 
recorded the instrumental settings directly. The operator 
applied the instrumental corrections to the cards by sorting 
the cards into groups and gang punching the appropriate 
correction. 

The conversion from corrected instrumental reading to 
star brightness was determined by an empirical process; 
this conversion was made previously by means of a calibra- 
tion curve for each plate. Each plate contained several hun- 
dred images whose brightness was to be determined and 
about forty comparison images of known brightness. Pre- 
viously, the known images were plotted and the results for 
the unknown ones read from the curve. The plotting and 
reading of these curves had been very laborious. For the 
punched card method investigation showed that, if a stand- 
ard calibration curve were used for all of the plates, the 
resulting errors for a given plate could be represented as ;a 
linear function of the abscissas. Therefore, it was decided 
to use a standard calibration curve in the form of a punched 
card file and to determine the two constants of the linear 
function for each plate. We determined these constants by 
the method of least squares, using an equation for each of 
the thirty or forty comparison stars. 

The remaining operations consisted of arithmetic opera- 
tions which were performed mostly in the parallel manner on 
simple equipment such as the 601 . The statistical analysis of 
the errors was primarily a matter of sorting and counting." 

As a final example of simplifying the punched card pro- 
cedure we shall describe a case where the longest way 
around is the shortest way home. Frequently a computing 
problem appears intricate and complicated, involving itera- 
tive processes and long sequential expansions. Yet, if the 
problem is turned around, the advantages of the punched 
card method can be applied directly, and the solution is ob- 
tained. The problem arises in connection with a method of 
radio navigation developed during the war and known as 
Loran. In this method each of two fixed stations broadcasts 
a signal which is received by the ship or plane. The receiv- 
ing instrument indicates the difference between the dis- 

aDetails of this work are given in Punched Card Methods in Sci- 
entific Computation, W. J. Eckert (Thomas J. Watson Astrtoomi- 
cal Computing Bureau, New York, 1940), Chapter X. 
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tances from the two stations; the problem for the navigator 
is to find the positions his ship could occupy in order to 
obtain the observed difference in distance. Since long dis- 
tances are involved, it is necessary to treat the earth as a 
spheroid, and when the equations are examined it is found 
that troublesome expansions are involved. However, when 
you turn the problem around and compute the distance 
from an assumed point on the spheroid to one of the broad- 
casting stations, the computation becomes simple and 
straightforward. The solution was to adopt a uniform grid 
of points covering the desired area and to compute by sim- 
ple parallel methods all the distances involved. The result- 
ing distances, recorded in cards, were used as a table, and 
the required positions were determined by inverse inter- 
polation. Although the grid consisted of a two-dimensional 
array of points, it was necessary to perform interpolation in 
only one coordinate, the other being held fixed. Since the 



data were in cards, it was possible to separate them into 
sections to facilitate the work. In some areas interpolation 
along a parallel of latitude gave greater determinateness, 
and in others interpolation along a meridian was preferable. 
The tabular interval was chosen to permit linear interpola- 
tion in most areas, and the critical areas could be removed 
for the use of second-order interpolation where necessary. 
Finally, the data were arranged with the sorter and listed 
in the most convenient order for the draftsman who was to 
plot them on a chart. 

In conclusion I wish to emphasize the value of examin- 
ing each problem as a job for the simplest equipment rather 
than for the most recent and most powerful equipment in 
operation. In the two and one-half years that the SSEC has 
been in operation we have been impressed by the number 
of problems proposed for solution on it that can be handled 
effectively on a standard accounting machine. 
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O N E of the most frequently performed groups of chem- 
ical engineering calculations is that involving fractional 
distillation columns. In a continuous distillation tower, a 
mixture of two or more volatile liquids is fed into the 
column. Depending on their relative volatilities, the compo- 
nents redistribute themselves on a series of trays or plates 
(sometimes these "plates" are merely theoretical equilib- 
rium points) above and below the feed. At each theoretical 
plate, equilibrium is reached between the descending liquid 
stream and the ascending vapor stream. At the top and the 
bottom of the column (and often at intermediate points) the 
partially separated portions are taken from the column. At 
the top, part of the vapor is condensed and returned as 
liquid reflux to the column; at the bottom, some of the liquid 
is boiled and goes up the column. 

The properties of each compound in the feed mixtures 
are known. The total composition is also known. The col- 
umn may operate at atmospheric pressure, or it may be 
working at high or low pressures. The chemical engineer 
desires to know the necessary amount of reflux, the height 
of the column (i.e., number of plates), the temperature and 
composition on each plate, and the composition of the top 
and bottom product. For a multicomponent mixture, these 
calculations can become quite complex and laborious. 

There is considerable variation in the detail used in per- 
forming these calculations. For example, if one wants only 
an estimate of the number of plates required for a particular 
separation, an empirical formula will provide this answer. 
It is also possible to perform an extremely rigorous calcu- 
lation, accepting no simplications or shortcuts. However, 
the most frequent choice lies somewhere in between. 

In our work we use the Lewis and Matheson method, 
which is quite well-known to chemical engineers. In this 
method, the composition of each plate is determined from 
the composition of the previous plate by combining the 
equations stating the equilibrium between the liquid and 
vapor leaving a plate with the material balance equations 
relating the compositions and quantities of liquor and vapor 



at any level in the column. This may be expressed as in the 
following four equations: 

1, For use above the feed in calculating up the column: 



P V 



D 



Xfi — Xfi — x p J Xq J 

2. For use below the feed in calculating up the column: 

p_v;_, w_ 

Xfi — Xfi—\ 'p T f ~r -^b T I 

3. For use above the feed in calculating down the col- 
umn: 

PL, DP 

4. For use below the feed in calculating down the col- 
umn: 



•^«, — -^1 



n+1 



PU^ 

p V 



Xh 



WP 
V p 



nome^ncIvATure; 

L — liquid flow leaving plate — above feed (mols). 
V — liquid flow leaving plate — below feed (mols). 

V = vapor leaving plate — ^above feed (mols). 

V = vapor leaving plate — below feed (mols). 
D = product leaving top of column (mols) . 

W = product leaving bottom of column (mols) . 
p = vapor pressure of a component (millimeters of 
mercury). 

P = total operating pressure (millimeters of mer- 
cury). 

Xn = the concentration of a component in the liquid 
leaving the wth plate numbering from the bottom 
(mol per cent). 
Xi,, Xq = concentration of bottom and overhead products, 
respectively (mol per cent). 
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SiMPUFYiNG Assumptions 

It is to be noted that the use of the vapor pressure of the 
pure component implies that the equilibria in the system 
follow Raoult's Law. This is a common assumption used in 
absence of data to the contrary; if more exact data were at 
hand, we believe it would be practical in many cases to cal- 
culate a suitable pseudo vapor pressure to use for a given 
system. It might be remarked at this point that the use of 
vapor pressure data, rather than alpha or k values, was 
chosen, because the same cards can then be used for a given 
compound in many different mixtures without further 
work. We also expect to use these same cards for other 
calculations such as flash vaporizations. 

It will be found necessary to make one further simplify- 
ing assumption in common use. This assumption is that the 
reflux ratios L/V are constant throughout the column. As 
will be brought out later, we expect to be able to handle 
variable L/V by expressing it as a function of temperature. 

Machine; Computation Me;thod 

In determining the composition of one plate from the pre- 
ceding plate, it is necessary to set up one of the above equa- 
tions for each component and solve it for the new concen- 
tration. Only when the set of vapor pressures is selected at 
the correct temperature does the sum of the calculated con- 
centrations equal 100 per cent. As a rule in calculating 
manually, an approximately correct set of vapor pressures 
is obtained by estimating the temperature of the plate. In 
that case, each calculated value is divided by the sum of the 
concentrations to yield a distribution totaling 100 per cent. 

In the machine method, advantage is taken of the high 
speed of calculation (four substitutions into one of the 
above equations requires eight seconds). A group of tem- 
peratures with one degree centigrade intervals is selected 
and a calculation made for each set of vapor pressures. For 
every set, the composition is totaled and the sum is in- 
spected for proximity to 100 per cent. The calculation yield- 
ing the closest set of values is taken as the correct one. This 
new composition is used to calculate the composition of the 
next plate and so on. Although many more calculations are 
made than are actually used, the rapidity of the machine 
more than compensates for the extra computation. 

While this operation has been performed using an IBM 
Type 602 Calculating Punch, there is no rea,son to believe 
that it cannot be performed with the IBM Type 602-A Cal- 
culating Punch or the IBM Type 604 Electronic Calculat- 
ing Punch. We have been able to calculate columns with up 
to four components on the 602 machine. Six component 
calculations could be calculated with the 602-A if all count- 
ers and storage units are provided. 

At first glance, the solution and testing of four sets of 
equations as above seems difficult for the machine to per- 
form on each card cycle. However, as several of the terms 



are constant for one portion of the column, they may be 
combined. Thus, the equations may be reduced to the forms: 

Xn = Xn-i • K ' p ± C 

for calculating up the column; 



'n+l 



-©-'•© 



for calculating down the column; 

where K = V/PL 

K' = PL/V 

^ D W 

C = XojorxiYT 

^, DP WP 

L = Xq -77- or jt'ft 



V 



V 



In a preliminary step, the constants are multiplied by the 
vapor pressures or their reciprocals as required. The equa- 
tions then reduce to the form: 

Xn = Xn+t • K" ± C 

where K" = K - p or K' {p\ 
C" = CovC^^ 



-G) 



The 602 may be wired to calculate four sets of these 
equations and to accumulate the sum of the results. The 
over-all procedure used in the calculation is determined 
when the chemical engineer fills out a form which contains 
sufficient information to enable the tabulating department 
operator to proceed with the calculation. The steps are as 
follows: 

Preparation of Vapor Pressure Cards 

These need be prepared only once for each compound 
used (unless systems are encountered requiring different 
pseudo vapor pressures). After calculation, the cards are 
filed for future use. As the cards may be readily duplicated, 
it is possible that prepared sets of vapor pressure data for 
common compounds will become available. The data are 
calculated by means of the Antoine equation: log^o^ = 
A — B/{T-]-C), although it is perhaps better to interpolate 
actual physical measurement of the vapor pressure. Briefly, 
the preparation of the cards begins with the gang punching 
of A and 5 on a reproduced set of {T-\-C) masters. These 
masters contain the reciprocals of T -f- 230 (the value of C 
used here) with the corresponding T (temperature) from 
-100 to 270°C. A-B [l/(r-H230)] is calculated quite 
simply. The results are merged with a set of logarithm 
cards, and both p and l/p are punched on each detail card. 
Thus, we have a set of 370 cards containing temperatures. 
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vapor pressures, and reciprocal vapor pressures. If several 
sets are prepared at once, the time required per set reduces 
to approximately twenty minutes. 

Preparation of Master System Cards 

A portion of the vapor pressure cards is selected which 
will safely include any temperatures encountered in the 
column. Cards are selected for each compound and are suc- 
cessively reproduced into a single set of system masters. 
These cards form a set of operating vapor pressures (and 
reciprocals) for the problem. They may be used for any 
problems involving the same mixture. Figure 1 contains a 
listing of the vapor pressure portion of cards for mixtures 
of carbon tetrachloride, trichloroethylene, 1,1,2 trichloro- 
ethane, and perchloroethylene. 

Preparation of Master Working Cards 

In this operation, the vapor pressures (or reciprocals) 
are multiplied by the appropriate constants. These cards 
are used for one problem only. 

Reproduce Detail Working Cards 

The master cards are reproduced into the cards that will 
be actually used in the plate calculation. Several sets of 
details are prepared. If additional cards are required, they 
may be reproduced later from the masters. 
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Figure 1 



Actual Plate-by-plate Calculation 

The working cards contain K", C", and the X7S that 
determines the sign of the C. Starting at the bottom, for 
example, the initial composition is gang punched into a 
set of details. The set is run through the machine, which 
stops and signals when the total composition passes 100%. 
The correct (closest to 100% total) card is selected and 
marked ;ri. This card is placed at the head of another group 
of detail cards, gang punched, and then set aside for the 
listing. This new group is calculated in the same way. This 
successive calculation continues until the terminating con- 
ditions are encountered. At this point, either a new section 
of the column is calculated or the calculation is discon- 
tinued. 



Printing the Results of the Calculation 

The correct cards are accumulated into a deck during 
the plate calculation with the bottom of the column on the 
bottom of the deck and the remaining cards in ascending 
order of the plate numbers. The cards are run through the 
405 accounting machine and listed. A partial listing for a 
distillation involving the same compounds shown in Fig- 
ure 1 is given in Figure 2. This listing is sent back to the 
chemical engineer as the results of the calculations. 

To illustrate the method using a well-known example, 
the three component system, benzene-toluene-xylene, was 
calculated using the conditions described by Robinson and 
Gilliland in Elements of Fractional Distillation. Figure 3 
shows the comparative results obtained by Robinson and 
Gilliland and by the use of this method. As these writers 
made exactly the same assumptions as are made here, the 
results should be comparable. While the feed plate, xylene 
disappearance plate, and total number of plates are the 
same, there are definite numerical discrepancies. The use of 
smaller temperature intervals (one degree steps instead of 
five) should slightly increase our accuracy. The three fac- 
tors discussed below, however, make the machine calcula- 
tions less accurate. 



Use of Calculated Vapor Pressure Data 

One cause of the discrepancy is the use of different 
sources of vapor pressure data. Robinson and Gilliland 
have read points from graphed physical data. Our vapor 
pressures were calculated by means of the Antoine equa- 
tion. 

In the range used for the calculation, the vapor pressure 
discrepancy had a mean value of 1.15% for benzene, 0.57% 
for toluene, and 2.25% for xylene. If the same data were 
used and correctly interpolated, this discrepancy would 
disappear. 
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Figure 3. Comparison oi^ Composition 
<y^ Be:nz:Ene;-ToIvUEne;-Xyi.e;ne; Syste;m 

(Robinson and Gilliland and by this Method) 



11 



11 



INDUSTRIAL COMPUTATION 



failure to Balance to Exactly lOO Mol Per Cent 

At each plate, the calculation is stopped at (usually) 
some value between 99 and 101 mol per cent. If each com- 
position were divided by the total, the distribution would 
remain the same, but the values used for the next calcula- 
tion would be more accurate. This division by the sum is 
perfectly feasible in machine operation but requires the use 
of a more elaborate calculating machine. 

Digital Errors 

The four multiplications are carried to six significant fig- 
ures, but two of the multiplicands are limited by machine 
capacity to the first three significant figures. This results 
in a digital error of .01% for two of the factors and .1% 
for two others. The use of a 602- A or another larger 
capacity machine would eliminate this inaccuracy. 

Because of the iterative nature of the calculation, the 
errors accumulate and are compounded. The nature and 
extent of these inaccuracies should be kept in mind. 

Once the required vapor pressure cards are on hand for 
the desired calculation, the time of calculation is as follows 
for a four-component calculation over a fifty-degree boiling 
range that requires forty plates for separation: 



Machine Operation 



Machine Time 



Prepare master system cards (1 set) 2.0 minutes 

Prepare master working cards (2 sets) 12.5 minutes 

Reproduce 5 sets of working cards (10 sets) 5.0 minutes 

Calculate plate composition 26.6 minutes 

Print column composition 0.5 minutes 



TOTAL 



46.6 minutes 



The number of cards per test will vary from ten to two. 
An average of five was used in calculating machine time. 

The additional operator manipulation time depends on 
the skill and experience of the operator. It should certainly 
not exceed the machine time. 

Calculations with Vanishing Components 

If one or more of the feed componeiits in the upper or 
lower part of the column is reduced to less than ,01 mol 
per cent, the calculation may be handled as follows: 

1. No vanishing components — Calculate up column or 
down column. 

2. One or more components vanish in upper section — 
Calculate up the column. 

3. One or more components vanish in the lower section 
— Calculate down the column. 

4. At least one component vanishes in each section of 
the column — Calculate both up and down the column. 

The practicability of applying machine methods of calcu- 
lation to chemical engineering calculation hinges on many 
factors. Rigorous calculations of this type should always be 



made by or under the direct observation of one who under- 
stands the problems involved. The method described here 
is offered as a rapid method of accumulating information 
about column design and performance. Some of the cases in 
which this rapid collection of information should be excep- 
tionally useful are suggested below: 

Calculation of columns for compounds and classes of 
compounds frequently handled, For organizations calcu- 
lating large numbers of columns, this should be quite 
valuable. 

Calculations in which several assumed conditions must 
be calculated in order to select the correct one. This is the 
case where the boiling point of one or more components 
lies between those of the components whose separation is 
desired. 

Preparation of correlations between properties, oper- 
ating conditions, and column performance. The effect of 
reflux ratio on column height, optimum feed condition, 
position for side streams, etc., may be found for a variety 
of systems. 

Calculation of temperatures and compositions for a 
number of variations in the operating conditions (e,g,, 
estimating optimum temperature control locations). 

In rigorous calculations the chemical engineer often 
makes a heat balance on each tray as he calculates down 
the column. This determines the change in L and V due to 
heat losses from the column, the change in liquid and vapor 
sensible heats with temperature, and the differences be- 
tween the heats of vaporization of the various components. 
Experience indicates that this change may be often closely 
approximated as a function of temperature alone for a given 
problem. This function may be found by first calculating 
the compositions for a constant L/V followed by thermal 
balancing at selected plates in the column. Interpolation at 
1°C. intervals will then give a first approximation to the 
L/V throughout the column. 

Since an individual detail card contains the vapor pres- 
sure and the reflux data for each temperature, it is clearly 
possible to use the interpolated L/F's on successive cards 
at different temperatures. By repeating the plate-by-plate 
calculations with the more exact L/F's, a more accurate 
calculation will be obtained. 

DISCUSSION 

Professor Donnell: The general advantage of doing the 
calculation by the punched card machine rather than by a 
slide rule is that very quick short-cut methods using 
abridged formulae yield insufficient information for most 
purposes. They may give one answer, the number of plates, 
but will not tell you temperature gradient or the composi- 
tion of a given point. 
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Professor Nichols: How do you go about calculating the 
feed plate in your position? 

Mr. Opler: We assume that we are given the composition 
that is entering the column at that point. Knowing that 
composition, we simply use it as our initial condition in a 
distillation. I actually insert a zero for the amount of the 
material which is removed at the top or bottom of the 
column. 

Pro jess or Donnell: I have a question in regard to the dis- 
appearance of compounds in the bottom and the top. How 
do you know the composition of the feed tray? 

Mr. Opler: We are assuming that the feed composition is 
the composition of the feed plate. This may be a weakness. 

Professor Donnell: Don't you think that is quite far from 
the case ? 

Mr. Opler: The trick is to match the feed as closely as 
possible to the composition of the feed plate. I am using a 
little inverse logic in trying to match the feed plate as 
closely as possible to the feed. 

Professor Donnell: The reason we are so particular is 
that we work in petroleum and have disappearing com- 
pounds in the bottom and the top. We have found that the 



composition of the feed tray is different; the reflux ratio is 
about four to one. You could go on and say, "We will just 
match the ratio of the key component on the feed tray with 
the feed. That does not give you the maximum feed tray 
location." 

Mr. Opler: That is very probably one of the weaknesses. 

The problem of choosing the optimum feed tray location 
is often solved when using the Lewis and Matheson method 
at the desk by trial and error, introducing the disappearing 
components as a chosen feed tray is approached from both 
top and bottom and carrying the calculation a few plates 
past the chosen feed tray. When the proper tray at which 
to introduce these disappearing components has been estab- 
lished by trial so that the compositions match sufficiently 
close at the feed tray, one completely balanced solution has 
been found. Repetition with choice of a different tray for 
the feed will disclose whether more or fewer total plates 
for the column are required. Further repetition will locate 
the optimum feed location for minimum total plates. Al- 
though we have not worked out examples of this in detail, 
we believe that the machine procedure would save consider- 
able time in carrying through the tedious repetitions. 
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READILY available large-scale computing machines 
such as the IBM Type 602-A Calculating Punch and Type 
604 Electronic Calculating Punch have been shown to be 
very advantageous for any calculation problem where multi- 
ple repetition of a relatively simple calculation procedure 
was necessary. For instance, tables of values of plate com- 
positions at finite reflux ratios for binary mixtures in 100- 
plate distillation columns, and for ternary mixtures in 
25-plate columns, have been readily obtained using the 604.* 
For more complex problems these machines are restricted 
by their limited storage capacity and by the fact that a new 
set of wiring for the control panels is necessary for each 
new calculation procedure used. 

The new IBM Card-Programmed Electronic Calculator, 
however, has approximately five times the storage capacity 
of the type 604, and a single set of control panels can be 
used for any arithmetical calculation procedure, provided 
fractional exponents are not involved. Thus, this machine 
overcomes many of the shortcomings of the previously 
mentioned machines. The card-programmed calculator owes 
its versatility to the fact that the control panels are wired 
to make possible each of the four basic operations of arith- 
metic. The particular calculation procedure desired is car- 
ried out by means of a deck of specially punched IBM 
cards, each card in its turn choosing the proper arithmetical 
operation and the proper factors to carry out one step in 
the calculation. A new equation can be evaluated merely 
by the use of a new set of program cards, properly punched 
to carry out the steps of the equation desired. Thus, by the 
proper combination of program cards any length problem 
can be handled, provided only that the storage capacity of 
the calculator is not exceeded during the calculation period. 

This versatility of the card-programmed calculator led 
the authors to investigate its use for the trial-and-error 
type of problem which arises in many chemical engineering 
design problems as well as in other engineering fields. 
Within the field of chemical engineering the unit operation 

* Solution of trial-and-error-type problems. 



of distillation was chosen because of its wide usage and be- 
cause of the large number and variety of problems en- 
countered. Also, in the field of distillation, as well as all the 
other diffusional operations, the trial-and-error problem is 
especially important, because in a large proportion of the 
cases studied, insufficient data are available to establish the 
value of all the independent variables present. Thus, the 
values of some quantities must be assumed, and the accu- 
racy of these assumptions must be checked — hence, the 
trial-and-error calculation. This situation arises almost con- 
stantly in the problems encountered in design studies on 
distillation columns. Solution of these problems by machine 
methods should result in larg6 savings in time and money 
for all organizations involved in this and related fields. The 
methods of solution of several variations of probkms of 
this type are discussed in this paper. 

ThiE Ge^ne^ral Probi.e:m 

A problem frequently arising in distillation involves de- 
termining either the head composition {xd), the number of 
plates (w), or the reflux ratio {R) corresponding to a par- 
ticular feed condition, with other column variables known 
or specified. In nearly every case the problem becomes a 
trial-and-error type of calculation, because an explicit 
function of either Xd, n, or R, in terms of the other vari- 
ables of distillation, is either very complicated or entirely 
inexpressible. This is especially true if the distillation sys- 
tem under study involves variation in relative volatility (a), 
heats of vaporization of the components, or plate efficiency. 

In order to define a distillation system properly, the val- 
ues of a certain minimum group of the variables must be 
established. This group need not always consist of the same 
variables, of course. The values of all other column vari- 
ables can then be calculated from these known quantities. 

For the case of an adiabatic continuous column (Fig- 
ure 1) operating on an ideal binary system (the simplest 
possible case), the minimum group of variables might con- 
sist of the following: 
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completely, can be calculated. The additional functions are: 

1. Enriching section liquid return rate (L) as mols per 
unit time. 

2. Stripping section vapor rate (F') as mols per unit 
time. 

3. Stripping section liquid rate (U) as mols per unit 
time. 

4. Bottoms takeoff rate (W) as mols per unit time. 

5. Bottoms composition (^Xw) as the mol fraction of the 
more volatile component. 

6. Feed plate composition (xi) as the mol fraction of 
the more volatile component. 

7. Plate compositions (jr„) as the mol fraction of the 
more volatile component. 

It must be kept in mind that variations of the problem in- 
volve other combinations of known and unknown factors. 

The usual design problem arises because not all of the 
minimum group of factors are known; and, in addition, the 
relations between variables are too complex to permit direct 
calculation. Therefore, the required procedure must be: 
(1) to guess at some probable values of the unknowns 
among the minimum group of factors; (2) use these trial 
values to calculate trial values of at least one of the de- 
pendent functions in two or more different ways; (3) from 
a comparison of these calculated values, determine a proba- 
ble error in the original trial values; and finally, (4) choose 
new trial values and repeat the above procedure until the 
chosen trial values show themselves to be correct. This trial- 
and-error procedure is important not only in distillation, 
but also in the design calculations of nearly every phase of 
chemical engineering. 



Figure 1. IdeaIvIZRd Diagram of a Typicai, 
Continuous Distii,i,ation Coi^umn 

1. Feed rate (F) in mols per unit time. 

2. Feed composition (sf) as the mol fraction of the more 
volatile component. 

3. Condition of feed (q) as a ratio of the heat required 
to vaporize one mol of the feed to that required to 
vaporize one mol of saturated liquid of the same 
composition. 

4. Distillate composition (xd) as the mol fraction of the 
more volatile component. 

5. Distillate takeoff rate (D) as mols per unit of time. 

6. Enriching section vapor rate ( V) as mols per unit of 
time. 

7. Relative volatility of the mixture (a). 

8. Number of theoretical plates («) in the column. 
From the above factors the following additional func- 
tions, which are necessary to define the distillation system 



Use of the Card-Programmed Calculator 
in Solving the Problem 

As a basis for a description of the use of the card-pro- 
grammed calculator in trial-and-error calculations, assume 
that the following problem is at hand: The variables whose 
values are known, or specified, are the same as those listed 
in the general case discussed above. Let the distillate com- 
position (xd) be the variable for which trial values are to 
be assumed until a satisfactory value is obtained. The com- 
plete calculation procedure is now carried out with the card- 
programmed calculator according to the steps of the follow- 
ing plan. The equations are numbered to correspond to the 
steps wherein they occur. 

1. The values of the known and assumed quantities are 
read into the calculator and stored in designated stor- 
age units. 

2. The bottoms rate is calculated as 



W = F - D . 



(2) 
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3. The bottoms composition is calculated as 
Pzf — Dxd 



Xw 



W 



(3) 



4. The enriching section liquid rate is calculated as 

L=V -D . (4) 

5. The stripping section liquid rate is calculated as 

U^L-^qP . (5) 

6. The stripping section vapor rate is calculated as 

V' = V+{q-\)F. (6) 

7. The feed plate composition is calculated as 



_ {L^-D)zf + D{q~\)xn 
(L+qD) 

8. The top plate liquid composition is calculated as 



Xt = 



Xd 



(X - {<x-l)xi, 



(7) 



(8) 



The immediately subsequent steps involve plate-to-plate 
calculations down the column by the well-known McCabe- 
Thiele procedure.^ That is, we calculate the vapor and 
liquid composition on each plate of the column in turn, until 
the bottom plate (reboiler) is reached. The liquid composi- 
tion in the reboiler is the bottoms composition, here desig- 
nated as x'w, if perfect mixing occurs. If the trial value of 
Xd were correct, then x'w would be identical with Xw from 
step 3. In general, the first trial value of Xd will not achieve 
this equality. A comparison of x'w with Xw gives a measure 
of the error in the original assumed value of Xd. From 
knowledge of this error we are in a position to choose a 
much more accurate value of Xd for the next trial. 

The plate-to-plate calculation of liquid and vapor com- 
positions by the McCabe-Thiele method is carried out on 
the card-programmed calculator as follows: 

9. The composition of vapor leaving the plate next to 
the top plate is calculated as 



L D 

yt-x =j^Xt-\-—XD . 



V 



V 



(9) 



10. The liquid composition on the plate next to the top 
plate is calculated as 



Xt-x = 



yt-i 



a - {ci-\)yt-i 



(10) 



Liquid and vapor compositions corresponding to the re- 
maining plates above and including the feed plate, that is, in 
the enriching section, are calculated similarly by repeated 
alternate use of the equations: 



11. 



L 



D 



yn — y ^n + 1 -\- y ^D 



(11) 



12. 



Xn — 



yn 



a - (a-1) yn 



(12) 



However, below the feed plate (that is, in the stripping 
section) a different set of equations must be used to cal- 
culate the values of yn and Xn because of the different flow 
characteristics of the column, caused by the introduction of 
the feed. The required equations are: 



and 



Xm, — 



V 

y, ^Wt + 1 — 


Wxw 

V 


ym 





(11a) 



(12a) 



a - (a-l)y^ 

The basis of determining whether a particular plate is in 
the enriching section or the stripping section of the column 
is to note whether or not Xn is greater or less than the 
value of Xi, the composition corresponding to the point of 
intersection of the operating lines with the q line. 

This can be used by the calculator to give it a basis of 
choice between equations 11 and 11a in calculating the 
value of y for the next lower plate. 

It can be seen that the equation 



Xn — {^i + 0.0001) = u 



(13) 



will give a negative value of u for any Xn less than or equal 
to Xi and a positive value of u for Xn greater than Xi. 

Use is now made of the negative balance feature of the 
machine to cause a selector to transfer when u becomes 
negative. This selector controls the region of the program 
card from which the program punches are read, and thus 
can cause the calculator to use equations 11 and 12 or 11a 
and 12a, if the proper control signals are punched in the 
two different regions on the card. 

The test described above must of necessity follow the cal- 
culation of every value of plate liquid composition and, 
therefore, will take place between calculation steps 10 and 
11 a.nd after step 12 as they are described above. 

To recapitulate regarding the plate-to-plate calculations, 
the calculator starts with steps 9 and 10 for the plate next 
to the top plate, then makes the comparison of equation 13. 
Unless the unusual condition of introduction of feed on the 
top plate is involved, the results of the comparison of equa- 
tion 13 will direct the calculator to proceed through steps 
11 and 12, following which the comparison of equation 13 
is again made. This cycle of steps 11 and 12 followed by 
the comparison of equation 13 is repeated until the test 
indicates the need for use of equations 11a and 12a instead 
of equations 11 and 12. These are then used repeatedly 
until the plate-to-plate calculations have been made for the 
specified total number of plates in the column. The value 
of Xm from the last use of equation 12a is, as mentioned 
before, the desired trial jrV- 
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When x'w has been calculated, the difference, 

xw - x'w = e-,^ , (14) 

is a measure of the error in the original assumed value of 
Xd. Because of the hyperbolic form of the equilibrium dia- 
gram for an ideal system, the function relating en^ and 
Cxj^ takes on the following form : 



C - 



(1.0-;r^)«=^. 



"w 



where 



K (x^r 



(15) 
(16) 



and K is a constant. 
From the value of C^j) calculated above, it is possible to 
determine a new value of Xd as: 

^D new trial value ^^ ^D i ^xj) ' \'-' ) 

Since the calculator must use integer powers of numbers 
in repetitive calculation, k must of necessity take the in- 
teger form, although this may require several more trial 
values of Xd because C^j) cannot be calculated exactly. 



Once the new trial value of Xd has been calculated, the 
whole procedure is repeated with the exception of step 1 to 
obtain a new trial value of Xjy, etc. The calculation is fin- 
ished whenever Ca,p equals zero, that is, the correct value of 
Xd has been assumed. Figure 2 shows two such trials as 
graphed on an equilibrium diagram. 

Other Variations o^ the; DistilIvATion Problem 
AND Their Solution 

Determining Reflux Ratios by Trial and Error 

As mentioned previously, the above-described problem is 
only one ramification of the general distillation problem. 
Another aspect which presents itself is as follows: A spe- 
cific distillate composition (xd) is desired, and the problem 
is to determine the proper reflux ratio (L/D) to attain 
this result from a known feed. This type of problem is at 
least as important as the one just discussed in detail and is 
solved in the same manner, except for two changes. In this 
case trial values of enriching section liquid rate (L) are 
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Figure 2. Graph to Show McCabe-Thiele, Trial-and-Error 
Method oe Solution of Distillation Problems 
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chosen, and a system must be devised to provide a basis 
for choosing new trial vahies of L. The incktsion of the 
Hquid rate (L) among the minimum group of required vari- 
ables, of course, places the vapor rate (V) among the de- 
pendent variables. 

The method used for the choice of the new trial value of 
L, and thus eventual determination of the proper value of R, 
takes the following form: 



R, 



l,/D =^ = U/Dr , 



(18) 



where r is > 1 when ^^^,, is positive and r is < 1 when 
ea.^y is negative. This condition is satisfied if 

r=l-i-K'{e,,,) , (19) 

and thus 



U = 



\+K'{e,,y) 



(20) 



Effects of Non-ideality 

Both of the problems considered previously have been 
for the case of a column in which all of the simplifying as- 
sumptions of distillation hold true. This was done in order 
to give maximum clarity to the trial-and-error aspects of 
the problem. A note is in order here, however, concerning 
some of the methods of handling the various types of non- 
ideality and to show that their only effect upon the problem 
is to lengthen it somewhat. 

The four most common causes of non-ideality in column 
operation are: 

1. Non-constant relative volatility (a) of the mixture 
being distilled. 

2. Non-adiabatic operation of the column itself due to 
heat gains or losses through the column walls. 

3. Unequal heats of vaporization of the components of 
the mixture being distilled. 

4. A column plate efficiency which is not equal to 100%. 
A discussion of the method of handling each of these will 
now be given. 

In most cases the components of the mixture do not ex- 
hibit a constant relative volatility over the complete com- 
position range; that is, the relative ease with which one 
component is boiled away from the other is not constant. 
This is due usually to physicochemical effects within the 
mixture itself. 

The method of solution involves one of two choices. 
Either the quantity a can be expressed as a function of x, 
usually a series, or x can be expressed directly as a function 
of y, again usually a series. Of these, the second choice is 
undoubtedly the best, for the other would involve an itera- 
tive calculation, because a is a function of x, and x is, of 
course, the quantity we desire to establish. In the second 
case, the calculations above would proceed the same as be- 
fore, except that equations 8, 10, 12 and 12a would be in 
the form 

•^'-/(3') , (21) 



different from and probably more complicated than that 
previously used. Equation 21 is commonly in the form of 
the series 



X = ay ± by"^ ± cy^ i 



(22) 



Non-adiabatic operation of the column has the effect of 
changing the values of L and V, and of course, U and V 
from plate to plate as one proceeds up the column. In the 
case of heat loss from the column, both L and V will be 
decreased by the amount 

Q 



and 



L -L -k-^ 



(23) 



(24) 



where Q is the amount of heat lost per unit time, and H^ 
is the heat of vaporization of one mol of the mixture being 
distilled. If the heat loss is constant for each plate, usually 
a good approximation, this merely means that steps 5 
through 8 must be repeated for each plate rather than occur 
only at the start of the calculation. In other words, the cal- 
culation for each plate now involves the use of equation 24, 
followed by steps 5 through 10 rather than merely steps 9 
and 10 or 11 and 12, 

The case of unequal heats of vaporization which also has 
the effect of varying L and V is best handled by the well- 
known Peters method^ which involves the use of a fictitious 
molecular weight for the components of the mixture so that 
H^ is again constant per mol of liquid. This has the effect 
of making the values of L and V, expressed in mols, con- 
stant although the values of L and V, expressed in pounds, 
may be decidedly different from plate to plate. The effect 
of the above change upon the calculations is usually only a 
change in the value of a. 

Theoretical plates — that is, column plates which fulfill 
equations 10 and 21 exactly — are uncommon in actual prac- 
tice. Therefore, an efficiency term must be inserted to com- 
pare the actual relation of x and y which exists in the 
column, to that of a column with theoretical plates. The 
efficiency term is usually expressed^ as 



B = 



^n+l 



^n + 1 



Xn 



(25) 



where .«•„* is the value of Xn expressed by equations 1 1 and 
12, and Xn is the value of Xn actually existing in the column. 
Therefore, 

Xn = Xn + i - B{Xn + l — ^n*) , (26) 

and the efficiency can be used in the calculations above by 
using equation 26 after step 12 for each plate. 

Consideration of the above discussion shows that for the 
general case any or all of the above causes of non-ideality 
can be taken into account in either of the calculations de- 
scribed previously with no major effect except to lengthen 
the calculations in some cases. 
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Gknerai, Sor^uTioN op TriaIv-and-Error Problems 

Consideration of the method of attack of the problems 
already discussed leads to some general statements regard- 
ing trial-and-error problems as a whole. 

The first and most important point to be brought out is 
that trial-and-error calculations can be done independently 
by the machine only if expressions such as equations 15 and 
20 can be derived. Other than this, the only restrictions are 
that arithmetical operations must be used, and the storage 
capacity cannot be exceeded. 

The general method of approach in calculating trial-and- 
error problems is as follows: 

1 . Express the known independent variables of the sys- 
tem along with a trial value of the variable to be 
established. 

2. Calculate the values of the dependent variables from 
known and assumed values of the independent vari- 
ables, being sure to express at least one of the de- 
pendent variables in two different ways, thus getting 
two possible answers for this variable. 

3. From the magnitude and sign of difference between 
the two possible answers of the dependent variable, 
determine a new trial value of the assumed independ- 
ent variable. 

4. Repeat the above steps until no difference is noted in 
the dependent variable as calculated by the two dif- 
ferent methods. The last trial value of the unknown 
independent variable is then the correct value. 

CoNCIvUSIONS 

The work of this paper has shown that the method of 
calculation control employed on the card-programmed cal- 
culator makes its use possible for nearly all types of trial- 
and-error calculations. The machine is capable of carrying 
out the calculations at a great saving in time compared with 
hand calculation, and this saving increases greatly as the 
complexity and length of the problem increases. 

The main requirements in applying the machine to 
common engineering problems where the trial-and-error 
method is involved are: 

1. The relations involved must be arithmetical or must 
be capable of being expressed arithmetically. 

2. A basis must be available to enable the machine to 
choose a new trial value of the desired quantity. 

3. The storage capacity of the machine cannot be ex- 
ceeded. 

Nomenclature; 

a,h,c, ... Constant coefficients used in equation 22. 

Ca;jy Correction to be applied to Xd to obtain a new 

Xd for the next trial. 
D Distillate takeoff rate in mols per unit time; as 

a subscript it refers to distillate. 
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B Plate efficiency as defined by equation 25. 

BiB^y Difference in the value of Xw as calculated by 

two different methods. 

P Feed rate in mols per unit time. 

H^ Heat of vaporization of the liquid mixture at the 

point in question. 
i As a subscript it refers to the intersection of 

operating lines, i.e., feed plate. 

K Constant in equation 15. 

K' Constant used in equation 20. 

k Exponent in equation 15, 

L Enriching section liquid return rate as mols per 

unit time. 
U Stripping section liquid rate as mols per unit 

time. 
m As a subscript it refers to the general value of 

the plate number in the stripping section. 
n As a subscript it refers to the general value of 

the plate number in the enriching section of the 

column, 

Q Amount of heat loss from column per unit time. 

q Condition of the feed as a ratio of the heat re- 
quired to vaporize one mol of the feed to that 

required to vaporize one mol of saturated liquid 

of the same composition. 

R Reflux ratio as L/D. 

r Divisor used in derivation of equation 20. 

t As a subscript it refers to the top plate. 

u Value of the difference between the liquid com- 
position on the plate in question and the liquid 

composition on the feed plate. 
V Enriching section vapor rate as mols per unit 

time. 
V Stripping section vapor rate as mols per unit 

time. 
W Bottoms takeoff rate as mols per unit time; as a 

subscript it refers to the bottoms. 
X Liquid composition as the mol fraction of the 

more volatile component; as a subscript it refers 

to the location in the column. 
Xn^ Value of the liquid composition as predicted by 

theory under ideal conditions. 
y Vapor composition as the mol fraction of the 

more volatile component; as a subscript it refers 

to the location in the column. 
Sf Feed composition as the mol fraction of the 

more volatile component. 

a Relative volatility of the mixture under study. 

A Symbol signifying the rate of change of the 

quantity in question. 
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Appendix 

Problem i. Calculation of the Correct Value of xj). 

I. The initial conditions are specified and supplied to 
computer: 

a = 2.23 F = 10.0000 Zf = 0.2000 

L = 4.0000 D = 1.0000 q = 1.0000 

11. The machine first uses equations (2) to (8) as in- 
dicated to calculate remaining initial conditions: 
W = 9.000 (2); xw = 0.1333 (3); F = 5.0000 
(4); X' = 14.000 (5); V = 5.0000 (6); Xi = 
0.2000 (7);xt = 0.2000 (8). 

III. The value of Xt is checked for feed plate location: 
M = 0.4420 (13) 

IV. Calculation of plates is then carried out in order as: 

(a) Plate 5: 

ys = 0.6737 (9.); X5 = 0.4808 (10); 
u = +0.2807 (13). 

(b) Plate 4: 

y4 = 0.5446 (11);. r4 = 0.3491 (12); 
w = +0.1490(13). 

(c) Plate 3: 

^3 = 0.4393 (11); .^3 = 0.2600 (12); 
u = +0.0599(13). 

(d) Plate 2: 

y2 ^ 0.3680 (11); .^2 = 0.2070 (12); 
u = +0.0069(13). 
(e) Plate 1: 

yi = 0.3256 (11); ^1 = 0.1780 (12); 
u = -0.0221 (13). 
The calculator would now begin to use stripping 
line equations (11a) and (12a), if there were more 
plates to be calculated. However, plate 1 is the last 
plate. 

V. Calculation of the new value of Xd is carried out 
using the equations (14), (15), (16) and (17) with 
Xd new trial value = 0.7666. 

VI. The above calculations are repeated for the new 
trial with the following results. 



b. Plate or 
Calculation 
Step 



Xw 



Jn 



Xn 



Xd 



Material Bal. 0.1378 

Top 0.5862 +0.3860 

5 0.6209 0.4234 +0.2233 

4 0.4906 0.3016 +0.1015 

3 0.3932 0.2251 +0.0250 

2 0.3320 0.1813 -0.0178 

Calculator transfers to stripping line equation here. 
1 0.2624 0.1376 

New;ri, -0.0002 0.7595 

0.7595 is, therefore, the true value of Xd- 



Problem 2. Calculation of the True Value of Xjj for a Case 
of Non-Constant Relative Volatility. 
Conditions are same as for problem 1, except that x is 
related to 3; by means of the equation below rather than 
eqtiation (12); 

X = ay + by"" = O.U29y + 0.8571/ (22a) 

Use of this relation also necessitated a change in equa- 
tions for the calculation of the new trial Xd so that 



■^D new trial value — 



xp + (LO-xpy 
SixwV 



(15b) 



This was necessary in order to damp the calculations 
properly because of the small spread in the equilibrium 
curve at its upper end (Figure 3). The actual operations 
are as follows: 

I. Initial conditions, both calculated and specified, are the 
same as for problem 1 except as specified above and 
below. 
II. The trials are carried out as follows: 



1. Plate or 
Calculation 
Step Xw 



yn 



Xn 



u 



-SBW 



Xd 



Material Bal. 0.1333 

Top 0.6628 +0.4627 

5 0.6902 0.5069 +0.3068 

4 0.5655 0.3549 +0.1548 

3 0.4439 0.2322 +0.0321 

2 0.3458 0.1519 -0.0482 

Calculator transfers to stripping line equation here. 

1 0.1854 0.0560 
NewxD +0.0773 0.8233 



a. Plate or 
Calculation 
Step 



Xw 



yn 



^xw 



Xd 



Material Bal. 0.1370 

Top 0.5956 +0.3954 

5 0.6298 0.4328 +0.2327 

4 0.4996 0.3093 +0.1092 

3 0.4008 0.2307 +0.0307 

2 0.3379 0.1862 -0.0139 

Calculator transfers to stripping line equation here. 

1 0.2748 0.1452 
NewjjTi) -0.0082 0.7595 



2. Plate or 
Calculation 
Step 



Xw 



yn 



Xn 



-WW 



Xd 



Material Bal. 0.1307 

Top 0.6985 +0.4984 

5 0.7235 0.5521 +0.3520 

4 0.6063 0.4017 +0.2016 

3 0.4860 0.2718 +0.0717 

2 0.3821 0.1797 -0.0204 

Calculator transfers to stripping line equation here. 

1 0.2679 0.0998 
NewA-i, +0.0309 0.8317 
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Equilibrium curve for 
the condition that . 

: 0.14293; + 0.857 ly 



x=y Line 




Equilibrium curve 
corresponding to 
a = 2.23 
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Mol Fraction of More Volatile Component in Liquid (x) 

Figure; 3 



0.9 



1.0 



3. Plate or 
Calculation 
Step ^w 



yn 



^^w 



Xd 



Material Bal. 0.1298 

Top 0.7117 +0.5116 

5 0.7357 0.5690 +0.3689 

4 0.6215 0.4199 +0.2198 

3 0.5023 0.2880 +0.0879 

2 0.3967 0.1916 -0.0085 



New 



1 

Xd 



u.jyo/ u.iyiD — U.UU55 
Calculator transfers to stripping line equation here. 
0.3028 0.1219 

+0.0079 0.8337 



4. Plate or 














Calculation 














Step 


Xw 


yn 


Xn 


u 


^x-w 


Xd 


Material Bal, 


0.1296 












Top 






0.7149 


+0.5148 






5 




0.7387 


0.5733 


+0.3732 






4 




0.6254 


0.4246 


+0.2245 






3 




0.5064 


0.2922 


+0.0921 






2 




0.4005 


0.1947 


-0.0054 







Calculator transfers to stripping line equation here. 
1 0.3119 0.1280 

New;fi> +0.0016 0.8342 



5. Plate or 
Calculation 
Step 



Xw 



yn 



^xyp 



Xd 



Material Bal. 
Top 
5 
4 
3 
2 

1 



0.1295 



0.7394 
0.6263 
0.5074 
0.4014 



0.7157 
0.5743 
0.4257 
0.2932 
0.1955 



+0.5156 
+0.3742 
+0.2256 
+0.0931 
+0.00^6 



New jTi) 



Calculator transfers to stripping line equation here, 
0.3143 0.1296 

+0.0001 0.8342 



0.8342 is, therefore the correct value of Xd. 



1. Warren L. McCabE and E. W. ThiELE, Industrial and Engi- 
neering Chemistry, 17, 605 ( 1925 ) . 

2. E, V. MurphrEE, Industrial and Engineering Chemistry, 17, 7A7, 
960 (1925), 

3. W. A. Peters, Industrial and Engineering Chemistry, 14, 476 
(1922). 

4. Arthur Rose and T. J. Williams, Industrial and Engineering 
Chemistry 42, 2494 (1950). 

DISCUSSION 

[There was an informal discussion of this paper during the demon- 
stration of the problem on the card-programmed calculator.] 



Application of Automatic Computing Nlethods 
to Infrared Spectroscopy 



GILBERT W. KING 
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ALTHOUGH it cannot be said that computing ma- 
chines have been indispensable to progress in infrared 
spectroscopy, they have allowed progress in certain fields, 
and as more people become active in the mixed field, ad- 
vances in almost all directions will be made. To illustrate 
both these points, we shall catalogue the instances where 
our computing laboratory has supported the infrared. 

Interpretation of Rotational Structure 

The earliest application was in the interpretation of spec- 
tra. The foundations for this are firmly based on quantum 
mechanics, and the dimensions of many symmetrical mole- 
cules have been deduced from the rotational structure of 
their infrared spectra. In the case of the majority of mole- 
cules, asymmetry prevented such an analysis because of the 
numerical work required. The most troublesome step is the 
calculation of the energy levels of a rotor with three differ- 
ent moments of inertia. These levels are simply related to 
the roots of certain matrices, and give rise to the well- 
known problem of finding characteristic values of sets of 
linear equations. In the case of the asymmetric rotor, the 
matrices (and their determinants) have only three non-zero 
diagonals, and the problem of finding the energy levels is 
simply that of finding the roots, rj, of a series of continued 
fractions. 



ao 



&i 



ai — 7] -\- &2 



= 0. 



(1) 



as — Tj . . . . 

This can only be done by successive substitutions of guesses 
at rj; however, a characteristic of continued fractions is that 
such a process rapidly converges. 

This problem is ideally suited to punched card equip- 
ment because the operation has to be repeated many times 
(spectroscopists would like a table consisting of a million 
T^'s) and because, further, the whole process depends on a 
simple arithmetic algorithm 



a — b/c , 



(2) 



which is iterated several times. These roots, incidentally, 
are the characteristic values of the Lame functions. 

A coarse table was constructed by these means. Five- 
point interpolation carried out on computing machines gave 
good approximations for rj at finer intervals. These approxi- 
mate ?;'s were substituted in the continued fraction to give 
the correct answers. In this way, a fundamental table of 
the energy levels of the asymmetric rotor has been built up. 
It has found great use by infrared and microwave spectro- 
scopists, and we are at present enlarging it. 

Having this basic table on cards, we can proceed with 
the analysis of rotational band spectra. The approach has 
been a stochastic one in which a structure of the molecule 
is assumed, the spectrum calculated and compared with the 
observed. Since several successive guesses have to be made, 
the repetitive nature of the procedure is suited to punched 
card techniques. 

The computation required is as follows: First, a para- 
meter, K, a characteristic of the molecule, is computed by 
hand, and the appropriate reduced energy levels, r]i{K), are 
looked up in the fundamental table. We then calculate 

Bi(K) = 7(/+l) [{a-c)rjM + (a+c)] , (3) 
where a and c are further parameters of the assumed dimen- 
sions of the molecule (actually the reciprocals of the least 
and greatest moments of inertia), and / is a second serial 
number labeling the levels. 

In the meantime, another basic table, giving the proba- 
bility % of a line occurring in the spectrum because of the 
absorption of energy transferring the molecule from energy 
level £{(k) to Bj{k), is taken for the k in hand, and repro- 
duced. This is a two-way table: First, the energies Bi{K) 
are collated with it, each Bi(K) card preceding all the ^^ 
cards with the same i. The value of Hi(K) is gang punched 
into the $y cards. Next, these cards are sorted on /, the 
Bj{k) cards are collated, and Bj{k) is gang punched. 
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We now have the basic cards representing the expected 
spectrum. On them we calculate the position of the line 

V = Bj(k) - EM , (4) 

and the intensities of the line 

ay = <I>,,-^-^i(«>/'^=^. (5) 

The exponential function of the argument is obtained by 
collating the $ cards with Hi{K)/kT = jit on them with a 
basic table of the exponential function e~'^. The latter is 
gang punched from the basic table into the cards. Multipli- 
cation by (^ij is an elementary operation. 

On sorting the cards by v, we have the basic line spec- 
trum, a card for each line, having its position v (in wave 
numbers) and intensity a. Unfortunately, as a rule, infra- 
red spectroscopy does not resolve the individual lines, and 
we come to a real contribution of punched cards. The re- 
corded spectrum is actually a certain integral of intensity 
over the finite slits of the spectroscope. We know the slit 
function p{v) well enough" and can calculate what our com- 
pletely resolved spectrum (the cards with a and v) would 
look like with finite slits by calculating the integral 



I{v) = 



J 



+ « 



p{v-v) (x{v')dv 



(6) 



The function of p{v) is replaced by a set of discrete values, 
and the integration obtained by summing over a finite in- 
terval. With card cycle transfer, this sliding integration can 
be done in one pass through the accounting machine. 

The above procedure is repeated many times for various 
structures until a satisfactory fit with the observed is ob- 
tained. The comparison is done graphically, and because 
many trials are made we have found our method^ of plotting 
with the 405 accounting machine very useful. Several infra- 
red bands have been analyzed in this way.^-^-^ Each one, a 
few years ago, could have been a doctoral thesis. 

Vibrational Spectra 

The majority of molecules have so many atoms that the 
rotational structure discussed above is smeared out and the 
infrared spectrum consists of many bands, each one asso- 
ciated with an excitation of vibrations in the molecules. An 
analysis of these spectra could proceed along the lines out- 
lined above, although the details are considerably different. 
Here, one would have to assume certain force constants in 
the molecule, calculate the vibrational frequencies by find- 
ing the roots of simultaneous equations, which in this case 
are more general than the ones discussed above, and recon- 
struct the spectrum from these results. Unfortunately, there 
are so many parameters to be introduced into the funda- 

«Recently we have carried out the necessary Fourier transforms 
giving the precise nature of p{.v), taking into account all the ex- 
perimental conditions, using punched card methods. 



mental equations, this procedure cannot be followed with- 
out further considerations. However, there is no question 
that computing machines will be necessary to solve the 
equations once a procedure for setting them up has been 
defined. 

At the present time, an empirical approach is being made 
in many laboratories to assign frequencies to certain struc- 
tural features. This is an empirical and, in a sense, a statis- 
tical procedure, and now calls for large-scale computing 
methods since several thousand spectra have been recorded. 
However, one should approach these data with considerable 
caution because there is no sense in trying to make empiri- 
cal rules out of incorrect data. In particular, the present-day 
infrared spectroscope is incapable of resolving all the fre- 
quencies which are characteristic of a molecule, and before 
any large-scale computing is done, the data should be better 
digested. One aspect of this that we are studying in our 
laboratory is the use of low temperatures to improve the 
resolution of the bands. Low temperatures decrease the 
natural widths of the bands and, therefore, separate them; 
but to see this decrease in width of the bands, one needs 
instruments of high resolving power. Computing methods 
to improve resolving power of an instrument are discussed 
below. To make the empirical assignments of band fre- 
quencies to structural features, the original data should be 
put on cards and comparisons made with the collator. The 
direct recording of spectra on cards is described below. 

Reflection Data 

Another example where present data are misleading and 
could give rise to a lot of fruitless calculation is shown by 
the reflexion spectrum of materials which are highly ab- 
sorbing, such as glass, quartz and other minerals. It is con- 
ventional to interpret the maxima in the reflexion curves as 
vibrational frequencies of the crystal. A closer examination^ 
of the theory shows that the reflexion curves are functions 
of both the absorption coefficient and the refractive index, 
and the characteristic vibrational frequencies of the crystal 
are determined from the maxima of the square of the re- 
fractive index times the absorption coefficient. It is, there- 
fore, necessary to get the refractive index and the absorp- 
tion coefficient from the reflectivity data. This can be done 
if the reflectivities are measured at two angles. In this case, 
the required parameters n and k can be obtained by a solu- 
tion of the transcendental equations: 
R = (7) 

(^ - 2f cos <^ + cos^</>) jg -f sin^<^ tan^c^) 

(c^ -f 2/ cos </) + cos2<^) (g + 2f smct> tan </> + sin^*^ tan^,^)' 

t>These experiments will be presented in detail elsewhere by H. O, 
McMahon and I. Simon. 
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1 



+^V[n^ (l-'«^) + sin2<^]2 + 4wV (8) 

[«2 {l-K^)-sm^] + V[«^ (l-K2)-sin2<^]2 + 4;iV 
and / = (9) 

Clearly, a solution for n and x, given two reflectivities* at 
two angles, is a considerably complicated problem. This has 
been done in our computing laboratory, and a fundamental 
table giving n and k in terms of two reflectivities has been 
made. Given two reflectivity curves, at each wave length 
one can enter the table by use of a collator and obtain the 
refractive index and the absorption coefficient. When these 
are plotted, it is seen that there is considerable difference 
between the true vibrational characteristics of the crystal 
and the observed reflectivities. Again, we note that the con- 
version of the recorded spectra into the final significant 
curves requires the reading of the signal at every point 
along the curve. Clearly, the direct recording of the data as 
digital numbers on cards, rather than as a graph, saves a 
great deal of manual labor. 



Card-Programmed Selective Sequence Calculations 

The solution of the equations given above involves a 
great number of operations and was only made possible by 
certain developments in computing methods with IBM 
equipment. Specifically, we have converted our 602-A cal- 
culating punch into a card-programmed selective sequence 
calculator. A single control panel is permanently wired. 
One of twenty-four different operations can be called for 
on this control panel by a lead card in which the required 
operation is coded. Other instructions have to be given to 
the machine, such as clearing counters and calling for a new 
card. The code numbers are punched on a lead card. Several 
successive operations can be calculated on each card by 
having a sequence of codes punched on the lead card. At 
the present time, we have only extended the calculations to 
three successive operations. With this control panel, we 
have found the 602-A more flexible than a desk calculator. 
With the use of this card-programmed control panel, the 
equations for the reflectivity are not formidable, and the 
number of operations required to evaluate them can be car- 
ried out in a reasonable time. 



INDUSTRIAL COMPUTATION 

Reading of the Data on Cards 

The foregoing examples are typical of the processing of 
experimental data by various theoretical and calibration 
procedures in which the recorded data have to be read 
point-by-point over a long strip. One spectrum alone may 
require 2,000 points to be read off in a continuous curve. 
Obviously, this is the bottleneck when a great deal of ex- 
perimental data are accumulated. To overcome this problem 
we have built an instrument, which instead of giving a 
continuous record of a signal taken experimentally, con- 
verts the signal into a number, which it punches in a card. 
Thus, the data are recorded as digital numbers on cards as 
well as on a continuous curve. Any processing of the data 
can be done with standard IBM equipment. For certain 
practical and theoretical reasons, the most efficient method 
is to record a single number in a single column of a card. 
To do this, the number has to be represented in a binary 
system. Fundamentally, punched cards allow a representa- 
tion of binary, rather than decimal numbers, for either there 
is or is not a hole in a certain location, corresponding to 
the digits 1 or 0, which is all that is required on the binary 
scale. We have found that the most convenient method of 
using standard IBM equipment is to indicate the binary 
digits by punches down the rows so that each column rep- 
resents a single reading of the instrument. A card then car- 
ries 80 successive readings. 

Our first application of this recording of an output of an 
infrared spectrometer as digital numbers on cards for proc- 
essing was a study of the noise of the instrument. This is a 
very appropriate example, because the noise is essentially 
discontinuous, and it is very difficult for a galvanometer 
continuous recorder to give an accurate account of noise 
arising in the detector. The digital reader does not suffer 
from the discontinuity of the signal. 

Contribution to Experimental Technique 

Our study of noise was the autocorrelation function to 
see if there are any fundamental periods which would be 
due to pickup, rather than to Jbhnson noise from the 
detector, 



<^(r) 



/T 



l(t-T) I(t)dt. 



(10) 



It is very easy to obtain autocorrelation functions on 
punched card equipment once the data are read off point- 
by-point as a digital number. One merely makes a new 
deck of cards, collates it into the primary deck, and sums 
all the cross products. The two decks are then separated, a 
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card removed from the second deck, and the two collated 
again. This is repeated, removing as many cards as re- 
quired, the sums of the cross products forming the auto- 
correlation function. 

To carry out calculations of this sort, it is not necessary 
to convert the binary numbers in the column of the card 
into a decimal number. By means of digit selectors on the 
602-A or the accounting machine, it is possible to read the 
binary number directly into the counters of the machine 
and carry out the multiplication on the decimal system. For 
example, if the value of the signal were 3 (in the binary 
notation 11), the cross product of 3 with itself would be 
read into the machine as 11 X H, giving 121 which is a 
mixed binary in decimal code, indicating 4 -\- (2x2) -{- 
1=9, the correct answer. 

The digital reader has also been used to record the spec- 
tra directly into cards. Autocorrelations of spectra have 
been made in order to determine the characteristics of these 
spectra for the design of filters; and cross-correlations of 
spectra with noise have also been made for use with a more 
advanced theory of filtering. Now that the data are in 
digital form and* can be handled by automatic computing 
equipment, a vast new field of processing data is available. 
For instance, it is now possible to apply any filter charac- 
teristics to the data without having to build special equip- 
ment. For example, it might be found desirable to use a 
square filter. The difficulty of building such a filter into the 
recording apparatus is that there are invariably phase shifts 
over such a filter which cause considerable unwanted oscil- 
lations of the recording equipment. When the values of the 
spectra are punched in cards, it is possible to process the 
data through the accounting machine and apply a filter of 
any characteristics without phase shift (or, strictly speak- 
ing, with a constant phase shift which is the same for all 
frequencies) . In particular, a square filter would be obtained 
by applying an equation like (6) with p{x) — sin x/x. 



We are at present developing the theory of time series to 
the recording of infrared spectra and anticipate employing 
various types of filters suitable for the spectrum in hand in 
order to improve the filtering and, therefore, improve the 
resolving power of the instrument. 

Analysis of Infrared Spectra by Use of Punched Cards 

It is quite clear when the experimental data are in 
punched cards in digital form, a great many other treat- 
ments of infrared spectra can be made. It is now possible 
to compare spectra not only at peaks, which has been the 
custom in the past, but in every wave length. One can also 
analyze spectra of mixed components by subtracting the 
desired amount of the spectrum of any one pure compound 
from the observed spectra. This would leave a residue 
which can then be studied for the presence of further 
components. 

Conclusion 

The application of computing machines to infrared spec- 
troscopy has been treated here very superficially, but it 
should be clear that punched card machines are now part 
of the scientific equipment of a research laboratory. 
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IT IS MOST unusual for any variable to be of great 
interest when considered by itself. Usually interrelation- 
ships of different quantities are sought in an attempt to 
explain what influence one might have on the other. The 
simplest type of interrelationship which can be imagined is 
one in which one quantity completely determines the other. 
Perhaps one of the most common illustrations of this rela- 
tionship is found in the arithmetic lesson of the grade school 
student who finds out the cost of three pencils if one pencil 
costs 5 cents. In a general formula form, if n is the number 
of pencils and c is the total cost, then the formula is 

c = 5n . 
Mathematicians call this a functional relationship, and it is 
said that c is a function of n because when n is given, c is 
completely determined. No doubt you can think of many 
others, especially in physics, because in this field many of 
the laws are expressed as precise functional relationships. 
It is for this reason that the use of statistics in physics is a 
rather new development. It has come into prominence in 
nuclear physics because of the possibility of electrons mov- 
ing in a variety of directions. It is now popular for the 
physicist to be interested in statistics. He must study some 
of the phenomena which are not simple functional relation- 
ships and which properly are described as falling in the 
field of correlation theory. 

Before attempting to define correlation, it might be well 
to consider the other extreme. The functional relationship 
is one of the boundaries of correlation, and this other ex- 
treme to which I refer is the other boundary. It can be 
illustrated by considering a somewhat senseless example. 
There is no relationship between the population of a state 
and the size of the shoe worn by the senator from that 
state. These two items are totally unrelated and they would 
be called uncorrelated. In the range from functional rela- 
tionship to totally unrelated variables there is a broad field 
which is included under the term correlation, although it 
must be remembered that it is common to think of the 
extremes as being special cases still within the frame- 
work of correlation itself. If the age of a married man is 
known, a reasonable estimate of the age of his wife can be 
made — although in most cases, probably, it cannot be con- 
firmed. This is due only to the reticence of most women to 



state their ages and not due to the lack of relationship be- 
tween the ages of husbands and wives. Except for recently 
notable examples which were in the newspapers, it is com- 
mon for husbands and wives to have approximately the 
same ages. It should be possible then to derive some for- 
mula for making a good guess, a guess which is based on 
the habits of the people and not purely on speculation. 

This illustration is exceedingly simple because it involves 
only two variables. This restriction is not necessary, and 
indeed it must be removed if certain problems in industrial 
statistics are to be solved. An excellent illustration of this 
is found in the doctor's thesis of one of my colleagues, in 
which he made a study of the demand for copper and in 
which an attempt was made to evaluate two methods of 
finding a relationship. In his analysis he tried to estimate 
the amount of copper delivered by looking at the undeflated 
price, the private gross capital formation, the stocks of 
copper at the beginning of the year, and the undeflated price 
of the previous year. This is a more practical problem and 
gives us a better illustration of the types of relationships 
which might be developed in actual practice. The important 
point is that these variables, as listed, are not sufficient to 
determine completely the deliveries of the copper, but they 
are highly influential and will include most of the items 
which combine to determine the actual deliveries which are 
made. 

In order to make these vague notions more precise, it is 
necessary to consider a few mathematical manipulations. 
It is best to consider first the simple case of only two vari- 
ables, because the fundamental concept for handling more 
variables is really no different from that involved in han- 
dling only two. Naturally, a search for the relationship be- 
tween two variables is started by looking at the values of 
these variables which have been observed in the past as 
occurring together. For example, find the ages of husbands 
and wives for 100 married couples and try to make an esti- 
mate of some formula from the data on the 100. Call the 
two variables X and Y, and assume that N is the number 
of pairs of values of these variables which are in the basic 
data. Consider X as being given and try to find a formula 
which would allow a prediction to be made of the most 
probable Y to associate with that X. Unfortunately, this 
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problem is very difficult if the formula is complicated. If, 
however, a simple formula is assumed, as in the case of the 
cost of the pencils, the problem is quite simple. Suppose 
that the relationship is 

Y = A^BX . 

The best such formula will be found, and what is meant by 
the word best should be defined. Obviously, it should be 
concerned with the size of the errors involved. By errors 
it is meant that, in an individual case, there is no guarantee 
that the Y value, estimated from the known X, is exactly 
correct. It is hoped that if there is an error in one case in 
one direction it will probably be counterbalanced later with 
an error in the other direction. This is not quite enough, 
since it is not desirable to allow a large error to creep in 
under the assumption that later there will be a large error 
in the other direction. The desired condition can be imposed 
by saying that the squares of the errors are to be small. 
This will force all errors to be small and eliminate the 
chance of counterbalancing large errors, one against the 
other. This criterion is called the criterion of least squares. 
If A and B are chosen by the principle of least squares, then 
the resulting formula is the one which will be called best. 
It might be well to keep in mind that if the relationship is 
at all valid and that X does form a good basis for the pre- 
diction of F by a formula of this type, this method or any 
other reliable method would give a close approximation to 
the correct answer. The pairs of values by X's and F's can 
be denoted with subscripts in the following way: 

{X^,Y,),{X^,Y^),....,{X^,Y^). 

For each X there is a corresponding estimate of F, and 
these X's and their estimates give N more pairs of values 
which will be written as: 

{x^,nx),{x^,n^),...,{x^,n^) . 

The following formula is assumed: 

F = ^ + 5X . 

An application of the least squares condition requires that 
A and B be selected so that 



I 



{Y.-BiY 



is less than it would be for any other values of A and B 
which might be selected. Substituting the formula for the 
B in this expression, it becomes 



I 



(F,- [A^BX,]y 



Take the derivative first with respect to the coefficient A, 
and secondly, with respect to the coefficient B, giving 



- 2y^{Yi-A-BX,) and - zV Xi{Yi-A-BXi) . 

i=l i=l 

The minimum value is obtained by setting each of these 
expressions equal to 0. If this is done and the summations 
are separated to simplify the form, the following result is 
obtained : 

NA + B^Xi - 2Fi,and 
A^Xi -f B^Xf = %XiY, . 

For convenience I have omitted the summation limits i= 1, 
i = N. Perhaps you remember your algebra sufficiently to 
see that it would be a simple matter to solve the two equa- 
tions for A and B with the result : 

, :^Yax\-:zxaXiY, 



B = 



N-^X\ - 

mXiY, 



■ ^xaYi 



mx\- (2X0= 



At this point, as is often the case in mathematical deriva- 
tions, it can be seen that it would be much easier if the sum 
of the Z's and the sum of the F's in these two expressions 
were each 0, for then A and B would have many fewer 
terms. Since the sum of the deviations of a variable about 
its mean is 0, this situation could have been obtained if the 
mean of X had been subtracted from each of the X values 
and the mean of F from each of the F values at the start. 
Denote these means by X and F, respectively, and let 

Xi = Xi — X and y^ = Yi — Y for all i's. 
Then "^Xi and ^yi would each equal 0. Furthermore, all of 
the steps above would be the same, but with small letters 
replacing capitals. The answers would have appeared in the 
simple form 

%Xiyi 



0, b 



%x\ 



this assumes that the original equation would have been 
written 

y = a -\-'bx . 

The best estimating equation would be 



y 



^Xjyi 



This is not a commonly used form. It shall be changed to 
the more usual one, and I will explain the advantages by 
the later discussion. This common form makes use of three 
constants: the two standard deviations a^ and oy and the 
correlation coefficient r^,, which can be defined as: 
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These can be introduced by the following manipulation: 

5^'f = Nal and ^Xiyi = Na^ayrasy . 

%Xiyi 



Therefore, 



N(Ta.CTyra.y 



O'v'^iC 



Substituting in the equation, we have 

CTyfxy 



y = 



X 



Of course, this equation can be expressed in terms of the 
original X and Y values by writing 



Y Y — ^'^'^^^ 



(X-X),or 



F= -^^^ iX-X) + F 



It was noted above that the errors in the estimates should 
be small. The size of these errors can be determined by 
calculating the standard deviation of the difference F — B. 
Denote this difference by d. 



Then d ^Y - U = Y 



Y -Y 



L 0-* 



(X-Z) + F 



1 



Ol/lX 



^^^ (X-X) =y- ---^x 



For convenience I shall omit summation subscripts, i, in 
the future. 



cri 



A^ 



2 2 

xy^- 



-H 



2^y^xyi-^xA 



X..2 O '^y'^^'BV AXy _. O-yKxy AX „ 



<Ty 



2 o ^y^xy i V (yx^y^xy i 



o-i/r. 



A^ 



2 



ary 2 



2 _ 9 2„ 2 I 2„ 2 
^y ^^y'xy ~r ^y'xy 



(T« — (T„r 



,2^ 2 



— 2/1 _ 2\ 

— CTy (^ i 'xyj • 

The standard deviation of these differences is called the 
standard error of estimate and is usually denoted by Sy. 
The formula for Sy shows that the constant r^.^ is never 



numerically larger than 1, since the quantity calculated is 
non-negative. The standard error of estimate is only in 
case the r^-y is equal to plus or minus 1. Before leaving this 
subject to discuss the more complicated cases, it is well to 
find an answer to the question, "What is the variability 
present in the estimates, B, themselves ?" This can be found 
by calculating the standard deviation of these B values. The 
mean of the B's, B, is: 



A^ L cr. 



(Z-Z) + F 



] 



N 
= Y . 



1 ^^^'^^2(Z-X)+^A^F 



N 



Then the variance of the £'s, al, is: 






al = ^%{B-By 






{X-X) ^Y -Y 



]' 



12 2 _ 



S;tr2 



A^ oi 


^r 2 1 

(Jy'xy i- 


ai N 


2^ 2 

Oyfxy 2 

;?— (T^ 



— ^y' xy • 

This says that oe = ^y'^xy • 

The standard deviation of the £'s, the standard error of 
estimate and the standard deviation of the original F values 
are all related by an interesting formula which can most 
easily be expressed in terms of the squares of these various 
standard deviations. Reference to the formulas just derived 
immediately shows that 

C2 2 2 

>^y = ^y ~ ^E , 

which is more usually written 

2 C2 I 2 

erg = Oj, -f- (T£! . 

This is the first time the statistics student encounters this 
formula, which is later termed the analysis of variance; 
the square of the standard deviation is frequently termed 
variance. In this case the variance of the original F values 
is broken into two parts, one of which refers to the vari- 
ability in the deviations between the estimates and the ac- 
tual values, Sy, and the other which refers to the deviations 
in the estimates themselves, o-f. It is common to speak of 
the latter as that amount of variability in the F variable, 
which is explained by the variable X. The former is that 
portion of the variability of F which is unexplained. It is 
clear that if Sy is large and o-f is small, the variability in F 
is not particularly associated with the variability in X. On 
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the other hand, if Sy is small and o-J is large, then a great 
deal of the variability in Y can be explained in terms of the 
independent variable X. The extreme case is that in which 
one of these two is 0. For example, if Sy is 0, then all the 
deviations are and every estimate is perfect. This situa- 
tion occurs in examples such as the one involving the cost 
of the pencils. On the other hand, if al; is 0, then none of 
the variability is explained in terms of X, and the situa- 
tion is of the type illustrated by the population of the 
state and the size of the shoe of the senator. Usually the 
intermediate situation arises and the attempt to explain Y 
in terms of X is partly successful. 

It might be well to consider a numerical example 
(Table I). For convenience, I have chosen one-digit val- 
ues for X and Y, and I have assumed AT to be 10 so that 
the division is easy. The arithmetic of the calculation is 
carried out in terms of the original variables, as well as 



in terms of the small letters which denote deviations from 
the mean. By the proper choice of values it was possible 
to obtain means which were whole numbers. This simpli- 
fied the arithmetic considerably. In this instance r^y was 
calculated to be .867, indicative of a relatively high degree 
of association between X and Y. 

It is not the purpose of this paper to discuss significance 
tests. It is easily seen that the importance of a result of 
this type is dependent on the number of pairs of values on 
which the calculations have been based. 

A very illuminating approach to the subject of correla- 
tion can be obtained by an examination of the relationship 
between the observed values of the dependent variable, Y, 
and the estimate of this variable, B. The correlation be- 
tween these two variables is exactly the same as the cor- 
relation between X and Y. This can be easily established 
if the small letters are used instead of the big letters. P 
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2 SO 40 314 212 250 

*There is no rounding in the tabular values. 

Z = 5, F = 4, AT = 10 

2 SZ2 



64 52 50 
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199.0625 
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N 
2XF 



- X^ = 31.4 - 25 = 6.4 
-Y^ = 21.2 - 16 = 5.2 



N 

XY^ 

N 



2.5298 



(Ty = 2.2804 



Z2F 



250- 5(40) 



N(Jx<yy 



10(2.5298) (2.2804) "57.69 



50 ^0.867 = ^ 



Normal Equations: NA -f B%X = %Y 10^ + 505 - 40 

A^X + 52X2 = ^XY SOA + 3145 = 250 
Solution F = 0.09375 -f 0.78125X 
199.0625 



(te 



rEY 
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2F£- 



- 16 = 3.90625 



<TE = 1.9764 



Y%B 199.0625 - 4(40) 



NapcxY 10(1.9764) (2.2804) 4.5070 



3.90625 „^- 

- = 0.867 = r^y 
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shall not bother to prove the following statement, but it is 
true that if each of the values of one of the variables in a 
correlation is multiplied by any constant other than 0, or 
if any constant is added to each of the values, or if both 
of these things are done, the correlation coefficient itself 
is unchanged. Since the difference between ^: and X is 
merely that X minus its mean is x, this will not change 
the correlation. A similar statement can be made with 
respect to Y. The truth of the above assertion about the 
correlation of Y with its estimate can now be established 
by calculating 



f'vM — fve — 






(TyCTe 



remembering the relationship between the estimate and x. 
This formula can be expressed in terms of x by substi- 
tuting 



e — ^ ^ X 



in the formula. This gives 



1 



Na-yO-g 
which can be reduced to 



^y'^^^x, 



'IIP, ' ! 






The only unknown term in this expression is the standard 
deviation of the estimated values which can be easily cal- 
culated as: 

1 2 2 

which reduces to 



2 r 

0-e = 



xy^y ^X 



N 



or 



2 _ V2 2 



If the above relation is substituted in the original expres- 
sion for the correlation, and is simplified by the following 
steps, the result is apparent: 



= r. 



'%xy 



%xy 



N<T:„r, 



x'xyf^y 



Nc 



= r. 



This formula states that the degree of agreement between 
X and y, as measured by correlation, is the same as the 
agreement between y and its estimate, when this estimate 
is based on x by the assumed formula. It is this approach 
which will be found most useful in the subject of multiple 
correlation to be discussed later. 

There are other properties of simple correlation which 
might be of considerable interest. One of these important 



properties is the fact that if y is estimated from x and then 
an attempt is made to estimate x from y, the two formulas 
are not identical. By this I mean that the equation 

Y = A-\-BX 

could not be solved for X to obtain the best estimating 
equation for guessing at an X which is associated with a 
given Y. Therefore, there are actually two equations: one 
minimizes the squares of the errors in one direction, and 
the other minimizes the squares of the errors in the 
other direction. These two equations will be the same 
only in case r is equal to plus or minus 1 . Geometrically, 
the graphs of these two equations are straight lines, called 
the regression lines. They intersect at a point whose X 
value is the mean of the X's and whose Y value is the 
mean of the F's. In the case of perfect correlation, r is 
numerically equal to 1; the two lines coincide, and all 
points lie exactly on the lines. In case of no correlation, 
the lines are perpendicular — one horizontal and the other 
vertical. 

I have spent a great deal of time considering the simple 
case involving two variables, because it is by analogy that 
this can be extended to multiple correlation in which any 
number of variables are being considered. It is conven- 
ient to think in geometric terms. The simple two- variable 
case, geometrically, is the problem of fitting a straight 
line to a set of points, which will all lie on the same 
straight line if r^y is equal to either plus or minus 1, and 
which will scatter widely about the line if r^-y is close to 0. 
The d values are the vertical distances from the points 
on the line to the observed values of F, and the E values, 
or the estimates of F, are the distances from the hori- 
zontal base line to the point on the line. 

If there are three variables, the geometric analogy is quite 
simple because it can still be thought of in terms of pictures 
which can be drawn. Unfortunately, it is difficult to visual- 
ize more than three dimensions, but the geometers who dis- 
cuss these subjects carry over the language of solid geom- 
etry into higher dimensions by speaking of hyperplanes 
instead of planes. For convenience the three-dimensional 
terminology will be used. The formula requires the estimate 
of one variable in terms of two others. It shall be assumed 
to be a first-degree expression similar to the above for two 
dimensions, such as 

Z == A -\- BX ^i- CY . 

Now Z is estimated from X and F; before F was estimated 
from X. Geometrically, this represents a plane in three di- 
mensions. Values of d and E may be discussed just as be- 
fore. The E value is the value of Z obtained by substituting 
a given pair of X and F values in the right-hand side of the 
equation. It is a vertical distance from the base plane to the 
approximating plane. The d value is the vertical distance 
from this point on the estimating plane to the actual ob- 
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served point. The only difference in the formulas between 
this and the previous case is that of the estimating equation, 

Z = A+BX + CY. 
The formula for d is identically the same, 

d = Z -E, 
with, of course, the substitution of Z in place of Y as the 
estimated variable. Of course, the formulas ior A, B and C 
will be derived by a slightly more complicated process, since 
there are more letters involved; therefore, there will be 
more equations. The criterion of least squares still gives the 
methods necessary to obtain these answers through the 
media of calculus and algebra. 

For the present those equations will not be discussed, and 
the correlation coefficient itself will be considered more di- 
rectly. Assume that an analogous standard error of estimate 
could be calculated. If the formula were completely worked 
out for the coefficients A, B and C, the given values of X 
and Y could be substituted and the estimates, E, calculated. 
Then each B could be subtracted from its corresponding Z 
and the values of the d's found. Then the calculation of the 
standard deviation of these d values would give the standard 
error of estimate. The important point is that the standard 
error of estimate does not involve any complication due to 
the existence of the two variables X and Y on the right side 
of the equation. It is still merely the standard deviation of 
a set of c^'s. Now this standard error of estimate could be 
assumed to be of a similar form: 

This R^ plays exactly the same role as r^ did in the previous 
case. I assume it is clear that I was forced to use z in the 
subscript in place of 3; as I did before, because z is now the 
variable being estimated. However, in the other instance I 
put two subscripts on the r. There is an interesting point in 
reference to the r, as well as a new one about R. Reference 
to the formula for r will disclose the fact that it is symmetri- 
cal with respect to X and Y. If every X were replaced by a 
Y and every Y were replaced by X in the formula for r, the 
answer would be unchanged. It has not been established, 
but it is a fact that with the three variables this symmetry 
is not maintained. We must differentiate in our minds be- 
tween the variable which is being estimated and the vari- 
ables on which the estimate is based. For this reason it is 
common to use a notation as follows: 

Rz.xY • 
You might naturally guess that the X and Y can be inter- 
changed in this formula without changing the value of the 
answer. However, you could not interchange the X and Z 
or the Y and Z without changing the value. 

The notion of multiple correlation, as incompletely de- 
scribed above, is the first and simplest of the extensions of 
the simple case. It is a measure of the association between 
a given variable Z and its estimate based on two other vari- 



ables X and Y. The inherent relationship between Z and X, 
also, might be desired, for example, if the effects of Y on 
each of the variables X and Z were removed. Stated in a 
more direct form it is: calculate the amount of X that can- 
not be explained in terms of Y and calculate the amount of 
Z that cannot be explained in terms of Y. Then the rela- 
tionship between these two amounts is the relationship 
measured by means of a partial correlation coefficient. In 
another form, it is asked, how much of the variability in Z, 
which cannot be explained in terms of Y, can be explained 
in terms of that part of variability in X, which also cannot 
be explained in terms of F ? When stated in this form, the 
formula and the arithmetic are almost obvious. What must 
be done is to find the differences between the values of Z 
and the estimated value of Z, if this estimate is based only 
on Y. The same must be done for X; namely, find the dif- 
ference in the actual value of X and the value that would be 
obtained by estimating X from Y alone. 

These two calculations will result in two variables which 
can be thought of in exactly the same way as X and Y in 
the case of simple correlation. The resulting correlation co- 
efficient is called the partial correlation between X and Z 
with the effect of Y eliminated or with Y being held con- 
stant. The latter description is probably the least desirable 
of the two. 

Now return to the problem of finding the estimating for- 
mula for a case with more than two variables. For conven- 
ience the formula 

Z = A+BX + CY 

referred to before, will be used. The extension to more 
variables will be very easy if this case is thoroughly under- 
stood. First express the deviations between the estimated 
values of Z and the observed values. Then square these de- 
viations, add them all and then choose the numerical values 
oi A, B and C so that this sum of squares of deviations is 
less than it would be for any other selection of the constants 
A, B and C. The expression for this sum of squares is 

2[Z- {A+BX-{-CY)Y. 
At this point, again, calculus must be used to obtain the 
resulting equations. The steps involved are merely partial 
differentiation with respect to the three coefficients A, B 
and C. These partial derivatives are as follows: 
-25 [Z- {A-\-BX^CY)] 
-2%X[Z~ (^+5X+CF)] 
-2SF[Z- {A+BX+CY)] . 
The minimum sum of squares is obtained when these deriv- 
atives are set equal to and the resulting equations solved 
for the unknown coefficients A, B and C. 

NA-^B-^X + C^Y = %Z 
A^X H- B%X^ + C^XY = ^XZ 
A%Y + 5SXF 4- C2F2 = %YZ . 
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We shall consider the formal solution a little later. Now I 
would like to point out the general formation of these equa- 
tions. Note that the first equation contains all the terms in 
the original predicting equation with a few summation signs 
and a coefficient N. The second equation is similar to the 
first one, except that each term involves an extra X, and the 
A'' does not appear. The third equation is similar to the 
second, except that it has an extra Y in each term instead 
of an extra X. If there were more than three variables, there 
would be more than three equations. The number of equa- 
tions is always equal to the number of variables. The addi- 
tional equations in the case of additional variables are 
formed like the last two mentioned above, in which an 
extra letter appears in each term. This means that there are, 
essentially, a set of equations with coefficients, which in- 
volve cross products of different variables and squares of 
variables all added together with some simple additions of 
the variables themselves. These simple additions always 
occur in the first equation and the first terms of the others. 
Again, notice that the whole set could have been materially 
simplified if x's, y's and ^'s had been selected instead of 
X's, Y's and Z's. This modification would have eliminated 
the first equation because all the terms except NA would 
have been 0. Thus, if the small letters are used, indicating 
deviations from the mean, for each variable, there would be 
only two equations when there are three variables. In gen- 
eral, the number of equations will be the same as the num- 
ber of dependent variables on which the prediction is based. 
In the present example this is two, x and y. There would be 
also only two unknowns, b and c (a would be zero). Notice 
that in terms of the little letters the equations are 

b'^x^ + c'^xy = Xxz 

b%xy -\- c%y^ = ^ys , 

and involve only sums of squares and sums of cross prod- 
ucts. These sums contain deviations from the mean, and 
hence are numerators in the formulas for correlation coeffi- 
cients and variances, the square of the standard deviation. 
For example, ^xy is the numerator in the formula for cor- 
relation coefficient 



selves would enter in the positions off the diagonal. This 
can be done in the following way: 



fxy — 



%xy 

N<Ta;(Ty 



For completeness, I shall also write the formula for a|, 

2 S-y 

In this formula the numerator occurs as one of the coeffi- 
cients in the two equations. It is, therefore, possible to ex- 
press all of the coefficients in these equations directly in 
terms of correlation coefficients and variances. The result 
of this method of expressing it is to obtain an array of co- 
efficients in which the variances occur in somewhat of a 
diagonal position, while the correlation coefficients them- 



5^ 

N 

:^xy 



^xy 
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23/2 



c = 



%xz 

N 



^ys 



N 

Tgiy(Tx(TyO -p (TyC = TygCyOg . 

This array of coefficients, when placed within the proper 
symbol bars, is called the variance-covariance matrix. The 
detail of the numerical solution of these equations will be 
discussed in another paper; so I shall not discuss that sec- 
tion of the theory, but say merely that the symmetry of the 
system of coefficients is the foundation for the development 
of a very simple solution routine, which is commonly known 
as the Doolittle method. The details of this method are un- 
important for the discussion, because what is of interest 
here is the form of the relationship which does exist. 

Indeed, further amplification of the theory of the normal 
equations needed to obtain the coefficients in the predicting 
equation is unnecessary, since the concept of multiple re- 
gression itself was attained without any necessity for finding 
the exact method by which the coefficients could be deter- 
mined. This point is particularly important if you are in- 
terested in the subject of graphic multiple correlation 
methods or in the subject of curvilinear correlation. The 
concept of the correlation coefficient as a portion of the 
formula for the standard error of estimate easily suggests a 
method for doing correlations graphically. An approximat- 
ing curve, not necessarily a line, can be drawn from an as- 
sumed equation or by a freehand fit to a given set of data. 
The deviations from this curve can be measured directly 
from the graph. The variance of these deviations will then 
be directly analogous to the standard error of estimate. 
The standard error of estimate can then be solved for the 
curvilinear correlation by means of the formula 

In this case r^ and s^ have meanings a little different from 
those in the straight line case. They are now measured from 
a curve. 



Regression in Terms of Independent Variables 

Consider a set of variables denoted by 

Xi, X2, Xs, and Z4 . 

For convenience the discussion will be restricted to the 
case of four variables; the extension to more is obvious. Let 
their deviations from their means be 

Xi, X2, xz, and X4, . 
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Let Xi be the dependent variable and consider the Hnear 
estimate of Xx from X2 alone by the least squares method. 
This can be written 

A'l = C12 X2 . 

The estimate of Xi from X2 is, in general, less reliable than 
the estimate from X2 and Xs. Since V2 and xs are usually 
correlated, use of only that portion of Xs not included by 
previously using X2 could be considered. This can be done 
by including, as an additional variable, the difference be- 
tween Xs and its estimate based on X2. Denote this by X3.2- 
The estimate is assumed to be the linear least square esti- 
mate. Then the estimating equation for Xi would be 

Xi = C12 X2 + C13.2 X2,,2 • 

This could be extended to include X4, in a similar manner, 
namely, by adding that portion of x^ not included in X2 
and Xs as estimated by a linear least square estimate. De- 
note this by X4.23. Then the estimating equation would be 

•*"l = ^12 -^2 + C13.2 ^3.2 + C14.23 -*'4.23 • 

This equation yields a set of normal equations with a single 
unknown in each equation. The variables on the right are 
orthogonal polynomials involving the original variables 
X2, Xs and X4,. The predictions resulting from this method 
are identical with those obtained from the usual regression 
method. Indeed, this procedure is the basis for an obvious 



method of graphic multiple correlation. It is of little value 
as a numerical tool because of the heavy arithmetic. It 
should be noted that the addition of other variables leaves 
previously calculated coefficients unchanged. 

DISCUSSION 

Mr. Bejareno: Could you give me a reference on correla- 
tion based on curve fitting which would be applicable to 
engineers' problems where they vary from a straight line, 
and it is apparent that there is a curve of some sort. The 
book known to me is Graphical and Mechanical Computa- 
tion by Joseph Lipka (New York: John Wiley, 1918). 

Dr. Sherman: In connection with curve fitting, I have a 
few remarks. One problem which comes up in engineering 
is that if you clock your data and try to lit it by means of a 
polynomial, there are cases in which the polynomial carried 
to rather a high degree is not a satisfactory fit, simply be- 
cause the data may, for example, fit a hyperbola. There is a 
good section in William Edmund Milne, Numerical Cal- 
culus (Princeton University Press, 1950) on this point. 
Another book which I have found useful in curve fitting is 
Deming's book. Adjustment of Bxperim,ental Data. General 
problems of curve fitting in which the variables x and 3; may 
be subject to error are answered, and good examples are 
given. 
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THIS CALCULATION, performed in the IBM 
Technical Computing Bureau, provided a complete numeri- 
cal analysis of the behavior of a pile when struck by a pile- 
driving hammer. The results indicated what stresses 
occurred within the pile and capblock, as well as the pene- 
tration of the pile into the earth. This is an example of the 
replacement of costly experimentation by economical cal- 
culation. 

Previously, using only key-driven calculators, it has 
been possible to study only partially a few isolated cases of 
the behavior of a pile under an impact. This was due to the 
high cost in time and money associated with each case 
studied. By employing the high speed and accuracy of IBM 
electronic calculating machines to perform this repeated 
formula evaluation, it is possible to study the behavior of 
numerous pile types. The cost of each case studied thereby 
will be substantially reduced. 

When piles are driven as a foundation for a building or 
other structure, the load that they will carry safely usually 
is determined by measuring the penetration under the last 
blow of the hammer and substituting this figure in a for- 
mula. Many different formulas are in use, and they vary 
widely in the answers they give. This calculation provides 
a means of mathematically testing the accuracy of these 
formulas when applied to various types of piles and various 
ground conditions. 

Furthermore, when the engineer has been asked what 
stresses the pile or the pile-driving equipment should be 
designed to withstand, he has often been at a loss for 
an answer. He has known only that certain practices and 
equipment have proved successful in the past. To do some- 
thing new has meant proceeding by the often costly process 
of trial and error. 

Because a pile is an object of considerable length, pile 
driving is a problem in longitudinal wave transmission and 
impact, the basic principles of which were first investigated 
100 years ago by the French mathematicians. Saint Venant 
and Boussinesq. These principles have been ably set forth 
by L. H. Donnell in ASME Transactions APM-52-14. 
However, the problem is a very complicated one. The 
method of solution offered herein is based on approximate 
integration using a step-by-step calculation, and is of gen- 



eral interest because it can be applied to other impact prob- 
lems. The calculation can be done by hand with a slide rule 
or desk calculating machine; however, modern electronic 
digital calculating machines are especially well adapted to 
make the calculation quickly and easily. 

General Case — Basic Theory 

For purposes of analysis, the hammer, pile, etc., will be 
represented by a series of concentrated weights separated 
by weightless springs. Subscripts m—\, m, m-\-\, etc., will 
be used to denote order of position, and superscripts «— 1, 
n, n-\-\, etc., will be used to denote order of time. 

Referring to Figure 1, let Wm-i, Wm and Wm+i repre- 
sent successive weights, and Sm.-i, Sm and ►S'm+i the springs 
underneath the respective weights. Let the initial positions 



etc. 



0„_, 



0„ 



0. + , 




Figure; 1 
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of weights Wm, etc., at the beginning of impact be indicated 
by 0,n, etc., and the initial lengths of the springs, by Lm, etc. 
Also, let 

£>w, etc. = instantaneous displacements (inches). 

ZvOT, etc. = instantaneous lengths (feet). 

Cw, etc. = 12[Zy^ — Lot] = Dm — Dm+i = instantaneous 
amount of spring compression (inches). 

Km, etc. = elastic constants for springs Sm, etc. = force re- 
quired to produce 1" of compression Cm, etc. 
(pounds per inch). 

Fm, etc. = instantaneous forces (pounds) resulting from 
'^m, etc. 

Vm, etc. = instantaneous velocities of Wm, etc. (fps). 

Rm, etc. = external forces, such as ground resistance, af- 
fecting the motion of Wm, etc. (pounds). 

Zm, etc. = instantaneous net force acting on Wm, etc. = 

Fm-l — Fm — Rm, CtC. (pOUnds). 

t = time (seconds). 

Also, let superscripts w— 1, n, n-\-\, etc., denote successive 
time intervals A^, so small that with negligible error it may 
be assumed that all forces and velocities remain constant 
during each time interval. Then, if by previous calculation 
or otherwise, the values of V and D are known for some 
particular time interval w — 1, all values for C, F, Z, V and D 
for the next time interval can be calculated as follows: 

Let F^~^ and D^'^, etc., represent definite known values 
for time interval «— 1, and let C^, F^, etc., represent defi- 
nite values to be calculated for time interval n. Then the 
following general formulas apply: 



rn _ nn—i _ nn— 1 

^m — -^m ^m+1 

pn T7~ fn 

Zn pn r>« r)» 
m — J^m — l "m -f^-ji 



F^ = V-m-^ + ^V 



n _ l/n-1 I ( 
m — 1^ m T^ I 



32.17AA 
Wm / " 



(1) 
(2) 
(3) 

(4) 



Dl = Dl-^ + ^D:^ = Dl-^ + (12A0 VI (5) 

AF above is evaluated by using the standard formula for 
change of velocity 

Force X Time g Ft 



{Vx-V2) 



Mass 



W ' 



where W is the weight, and g is the acceleration of gravity. 
Al> above is evaluated by the formula for distance traveled 
s = vt with the coefficient 12 introduced to convert to 
inches. 

Careful consideration of the above formulas will disclose 
that the force at the beginning of an interval is used to cal- 
culate the velocity at the end of the interval, and then this 



end velocity is used to calculate the distance traveled dur- 
ing the same interval; therefore, the forces, velocities, and 
displacements used are slightly out of step with one an- 
other, depending on the size of the time interval. 

These five formulas are used repeatedly until the calcu- 
lation has been carried as far as necessary. They apply 
whether or not the successive weights, elastic constants, 
and external forces are equal or unequal. Furthermore, the 
values of K, R, and A^ may be changed from interval to 
interval according to any definite formula, or suddenly as 
required. Sudden changes are required, for instance, when 
the stress in a material reaches the yield point, when a 
weight loses contact with a spring designed only for com- 
pression, when a coefficient of restitution is introduced, or 
when a ground resistance force must be made negative so 
as to resist temporary upward movements. Such conditions 
are called boundary conditions. Some of the more compli- 
cated digital calculators can take care of these boundary 
conditions automatically. 



Choice of Lengths L and Time Intervals M 

It should be borne in mind that formulas (1) to (5) in- 
volve the use of small but finite increments, not infinitesi- 
mals. The lengths L and the time intervals Af must, there- 
fore, be chosen small enough to suit each particular type of 
problem. For each new type of problem halving or quarter- 
ing the size of the units and recalculating part of the prob- 
lem must be tried until it is found that the use of smaller 
units makes a negligible difference in the peak stresses that 
occur soon after impact. The time interval can be changed 
while a calculation is in progress by inserting a new value 
for A^ in formulas (4) and (5), although this change is 
subject to the limitation pointed out in Step 2 below. 



Illustrative Problem 

A pile of non-uniform section as shown in Figure 2 is to 
be driven through water or very soft mud to a hard layer 
of ground which is capable of resisting a maximum force of 
600,000 pounds under the pile point. If the point of the pile 
starts to move upward momentarily, a negative frictional 
force of 100,000 pounds must be assumed, acting to hold 
the point down. No other side frictional forces are to be 
considered. The calculation is made to determine the final 
penetration per blow at which the assumed point resistance 
of 600,000 pounds will be developed. 

Side friction along the pile has been omitted from this 
problem so as to allow the stress wave to travel with the 
known speed of stress (or sound) in the pile material and 
thus provide one way of checking the calculation. The 
method would apply equally well no matter what values 
were assigned to side friction. 
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Step I. Decide on the time interval A^. From previous 
experience a time interval of 1/4000 second has been chosen 
as being small enough to give accuracy within about S%. 

Step 2. Decide oh lengths L. These must be at least as 
great as the distance stress will travel in the chosen time 
interval A^; otherwise the stress wave will run ahead of the 
calculation, and the results will be meaningless. It is recom- 
mended that L be made equal to twice the distance that 
stress would travel in the chosen time interval. The upper 
part of this pile is entirely of steel, and the known speed of 
stress in steel is 16,800 fps; therefore, the recommended 
length for L is (16,800 X 2) -^ 4000 = 8.4 feet. The pile 
length, plus a little added for the follower, happens to be a 
multiple of this figure; therefore, 8.4 feet can be used 
throughout. If an odd length were required, it would be in- 
serted at the point of the pile. 

Step J. Prepare a diagram as per Figure 3 showing how 
the ram, capblock, follower, and pile are to be represented 



for purposes of calculation. The individual weights Wx, 
etc., are calculated so as to give a weight distribution 
closely equivalent to that of Figure 2. 

Step 4. Prepare a tabulation of all constants required for 
formulas ( 1 ) to ( 5 ) as per Figure 4. This is readily done 
by considering the weight, cross-section, and modulus of 
elasticity of each portion of Figure 2 equivalent to a single 
spring or weight in Figure 3. The elastic constant K^ for 
the wooden capblock must be determined by experiment or 
must be assumed. For this problem, a value of 6,400,000 
pounds per inch has been assumed, which represents a 
rather stiff capblock, perfectly elastic. 

Step 5. If the work is to be done by hand, it is con- 
veniently tabulated as shown in Figure 5, which covers 
only the first three time intervals. In order to start the 
calculation, it is necessary to have a value for V and a 
value for D in the first time interval. The value of V, 
representing the velocity of the ram at the beginning of 
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(At = 1/4000 Second) 
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7,206,000 
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695 
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14 
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See Not? 
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Note: R^ assumed=F,, until F^^ reaches 600,000§. Thereafter Ru=600,000f if 
Wi4 is moving down and— 100, 000 jl^ if Wn is moving up. 

Figure 4 
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the impact, must be calculated by considering the dis- 
tance it falls and allowing for hammer efficiency. For this 
problem, efficiency was assumed to be 90%, which gave 
a velocity of 14.25 fps to be used in starting the step-by- 
step calculation. The numerical value of the displacement 
D in the first time interval is obtained from the assump- 
tion that the ram continues to move with undiminished 
velocity through the first 1/4000 second after the impact. 
For an impact velocity of 14.25 fps this gives a displace- 
ment in the first time interval of 0.04275 inches. This 



displacement then is used to calculate force Fi for the 
second time interval, and so on. As the stress wave 
travels down the pile, additional columns are needed in 
groups of five, all headed by the basic formulas (1) to 
(5) using constants taken from Figure 4. 

If the work is done by an electronic digital calculator, 
the results will be tabulated by the machine as shown in 
Figure 6. This is a condensed tabulation which shows 
time intervals 1 to 5 and some of the later time intervals. 
For the later time intervals, only the data for the top 
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Subscript 
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Forces 
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Forces 
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0.4408 
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-0.5238 

0.1696 


0.66504 
0.66926 
0.30196 
0.28940 


59 


1 

2 

13 

14 


-27,017 
15,191 
246,113 









600,000 


-7.4977 
-2.5485 
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and bottom of the pile are included, because these are 
the controlling factors in the calculation. 

The calculation begins in the first time interval with 
a ram velocity of 14.25 fps and a displacement of 0.04275 
inches. The stress wave travels down the pile, bringing each 
successive section of the pile into action. For example, in 
the fourth time interval, the velocity F3 amounts to only 
1.5946 fps with a displacement of D^ of only 0.00591 inches, 
and a total force F^ of only 8100 pounds. 

It will be observed that, to represent the specified ground 
resistance, i?i4 is given a value equal to i^i3 until h\z ex- 
ceeds 600,000 pounds, as it does first at time interval 30. 
From this point on, i?i4 remains at 600,000 pounds, except 
that when F14 becomes negative, as it does in time interval 



56, i?i4 is given a value of —100,000 pounds in the next 
time interval, which, in this case, is interval 57. In interval 

57, the velocity F14 is again positive; therefore, i?i4 is again 
given a .value of 600,000 pounds in interval 58, and so on. 
In the meantime, the values of the displacements Dn have 
increased gradually to a maximum of 0.28940 in interval 58, 
after which they remain practically unchanged. It will also 
be observed that the capblock force F-^ becomes negative in 
time interval 59. This means that the hammer ram has 
separated from the capblock; therefore, at this point, the 
hammer ram is dropped from the calculation. 

Figure 7 is a graphic representation of the results of the 
calculation. A separate curve is plotted for each weight W , 
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and the curves represent the relationship between displace- 
ment and time. It can be seen from the curve for Wi that 
the hammer ram curve crosses the pile head curve in the 
58th time interval. On the curve for W14, it can be seen that 
the penetration reaches a maximum in the 53rd time inter- 
val and then fluctuates slightly in the succeeding time in- 
tervals, reaching practically the same maximum again in 
interval 58. The calculation should always be carried beyond 
the first maximum of penetration in order to make sure that 
only slight fluctuations in penetration will occur thereafter. 

Checking the Calculations 

The total energy of the system for any particular time 
interval can be obtained by adding the kinetic energies of 
the individual weights, the potential energies of the indi- 
vidual springs, and the total work performed in overcoming 
the various external forces. The total should equal the 
energy of the ram just before impact. 

The total momentum of the system for any particular 
time interval can be obtained by adding the products of 
each mass multiplied by its instantaneous velocity and the 
products of each external force multiplied by the total time 
it has acted. The total should equal the momentum of the 
ram just before impact. 

Neither check is exact because of minor inaccuracies 
in this method, but a sudden variation between one time 
interval and the next indicates a numerical error. If the 
total varies by more than about 5%, consideration should 
be given to reducing the time interval and possibly the 
lengths L. If the work is done by hand, it is recommended 
that checks be made at every 10th interval. If the work is 
done by an automatic calculator, it may be possible to in- 
clude a running check as part of the setup. The energy 
check is to be preferred to the momentum check as it is 
more complete. Plotting the calculated results for displace- 
ments D as per Figure 7 is also an excellent check on the 
reasonableness of the results. Curves for separate calcula- 
tions may be readily compared. 

Recalculation for Change in Ground Resistance 

A change in the resistance near the point of the pile will 
change only a half or a third of the total calculation. Piles 
may, therefore, be recalculated for various point resistances 
with a considerable saving of effort as compared with aii 
entirely new calculation. 

DISCUSSION 

Mr. Sheldon: I would like to make a few comments. I 
have been much impressed with Mr. Smith's handling and 
understanding of the phenomena which go into these piles. 
He has not employed calculus, but he has really derived for 
you the partial differential equation for the motion of the 
pile with boundary conditions. 



The equation for the displacement D(x,t) in a one- 
dimensional, inhomogeneous elastic medium is: 

where p{x) is the mass density, Y(x) is Young's modulus, 
and f(x,t) is the applied stress. The simplest difference 
system by which we can replace the above differential equa- 
tion is, in the notation of Mr. Smith, 



Pm 



(Aty 

ym + l[Dm+l ~ Dm] ~ Ym\Dm ~ D m—l\ 



+ fm. 



The solution of this difference system is exactly equivalent 
to the solution of the difference system derived by Mr. 
Smith from first principals. To see this, note that 



Km 






w 

32.17 



= Pm^x and Vm — 






\2m 



Mr. Smith chose an interval A^ = A.r/2c (c = sound 
speed), so that he had a safety factor of 2 in the Courant 
condition for the stability of numerical integration of hyper- 
bolic type equations. 

We solved this problem on the card-programmed cal- 
culator at the technical computing bureau. We chose the 
CPC for solution because Mr. Smith had quite complicated 
boundary conditions imposed. For one thing, the resistance 
of the ground is a non-linear function. It is 600,000 pounds 
upwards when the last weight is moving down and 100,000 
pounds downward when the last weight is moving up. An- 
other condition that has to be provided for is that the cap- 
block and follower are not attached to the pile itself, so 
that after a certain period of time the capblock and follower 
fly off the top of the pile. Using the card-programmed 
calculator with its facility to list answers as we go along, 
we were able to observe the sign of the velocity at the 
bottom of the pile and also whether there was a tension or 
compression in the capblock. As soon as the condition 
which bounded this motion changed, we were able to insert 
a new instruction card which would take care of the new 
condition. This is a much simpler procedure than attempt- 
ing to put these conditions into the control panel of a 
machine. We have solved, totally, eight cases of this pile- 
driving work, and for the last six cases there was an addi- 
tional complication in the auxiliary conditions. Mr. Smith 
decided to take account of the fact that the capblock was 
made of wood and, therefore, was not perfectly elastic. We 
changed the elastic constant of the capblock according to 
whether the capblock was being compressed or was ex- 
panding. 
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The problem runs at about one step in the time every 
three minutes. This is an avera:ge figure, taking account of 
changing the program cards to take care of auxihary con- 
ditions. 

Dr. Aronofsky: I do not understand the boundary condi- 
tions at the top. The weight of the ram is 7,500 pounds. Is 
there any condition imposed? 

Mr. Sheldon: The weight at the top is a freely falling 
weight; so the boundary condition at the upper end of the 
pile is that at ? = the ram has a certain definite velocity, 
and all the other weights are not moving. 

Dr. Aronofsky : Is there any assumption about resistance 
along the lateral side of the pile all the way down ? 

Mr. Sheldon: In one case there was just the resistance at 
the bottom. In another case the resistance was applied in 
the middle. 

Mr. Smith: In Figure 4 the only resistance that is in- 
serted is the last one. All the other resistances are zero. 



However, if desired, you can put in as many resistances as 
there are weights. 

Dr. Buchhols: I think such studies have been made on 
analog equipment. You replace this type of system by a 
network of little capacities and provide certain nonlinear 
elements to take care of boundary conditions and special 
conditions. You run into a bit of a problem in the case of 
the capblock leaving the rest of the system. I don't know 
whether this problem has been done, but I imagine it might 
be possible to do so. 

Mr. Moncreiff: Was special wiring used for this problem ? 

Mr. Sheldon: No effort was made to change the standard 
setup at all. We made it very simple so that it took about 
one day to plan for the machine. 

Mr. Moncreiff: The calculation is simple enough so that 
you could save time by wiring a special control panel. 

Mr. Sheldon: That is true. 
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THEEFEICIENCY with which punched card tables 
can be constructed and used is so powerful that each prob- 
lem should be examined for the extent to which tables can 
be used. Examples of tables commonly used in this manner 
are trigonometric and logarithmic tables. The cards are an 
ideal medium for mathematical tables; it is possible to 
look up millions of digits in a single day. In fact, the effi- 
ciency of the method increases as the number of cards 
increases. 

An ordinary printed table has a tabular interval, argu- 
ments, function values, and it may have interpolating aids 
such as differences. The simplest table is a critical or turn- 
ing point table in which the arguments are so chosen that 
the value of the tabular function changes by one unit. For 
example, consider a section of the sine table: 



TABI.E I 



Sine;r 



2°.551 
2°. 608 
2°. 666 
2°. 723 
2°. 780 



0.045 
0.046 
0.047 
0.048 



For values of the argument between 2°. 608 and 2°. 666, the 
value of the sine is 0.046. If the required value of the argu- 
ment may be one of the printed values, then the author must 
state which of the two adjacent function values corresponds 
to that argument. This is usually stated near the printed 
table. In this example, the lower of the two adjacent values 
is taken, that is, sin 2°. 608 is 0.046. This critical table of 
the sine consists of 1001 cards for the first quadrant, and 
the error is 0.5 in the third decimal. 

In most cases the number of necessary entries in a critical 
table is prohibitive; we may shorten the table by making use 



of the differences of the function. Consider the following 
small section of the sine table: 



Tabi.^ II 
Sine Xi 



8'i 



2°.0 


0.035 


+ 17 


3°.0 


0.052 


+ 18 


4°.0 


0.070 


+ 17 



5°.0 



0.087 



To compute the value of sin 3°. 38 using linear interpola- 
tion, we have 



/« = /{ + w 8i where u 



-^i+l ^i 



sin 3°.38 = 0.052 -{- ^ (0.018) = 0.059 . 

The required argument serves two purposes in a table in- 
volving interpolation; it is used to enter the table, and it is 
used in the interpolating process. This sine table involves a 
linear interpolation, consists of 91 cards, and has the same 
accuracy as that of the critical table. 

The number of entries in the table may be shortened still 
further by the introduction of higher order differences 
(Table III). 



Xi sm Xi 



81 



'i+l/2 



Tabl^ III 

8'i ( Si— i/2+8i+i/2)/2 8'{/2 



0° 


0.0000 


+ 1736 





1736 





10° 


0.1736 


1684 


- 52 


1710 


- 26 


20° 


0.3420 


1580 


-104 


1632 


- 52 


30° 


0.5000 


1428 


-152 


1504 


- 76 


40° 


0.6428 


1232 


-196 


1330 


- 98 
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Tabliv III (Continued) 
Xi smxi 8i+i/2 §7 ( 8i_ 1/2+ 8i4. 1/2)72 §7/2 



50° 


0.7660 


1000 


-232 


1116 


-116 


60° 


0.8660 


7Z7 


-263 


868 


-132 


70° 


0.9397 


451 


-286 


594 


-143 


80° 


0.9848 


152 


-299 


302 


-150 


90° 


1.0000 




-304 





-152 



Using Stirling's interpolation formula, we have: 



+ «kiz^ + ^|.1 



To compute sin 37°. 689, 
sin 37°. 689 = 0.5000 

+ .7689 p^+'^^« + 7689( - :°^)] = .6111 , 

The necessary interpolation procedure involves three addi- 
tions, two multiplications, and two divisions. In many prob- 
lems, however, it is more efficient to tabulate quantities 
other than the differences themselves. If the quantities 
( 8/- 1/2 +Si4. 1/2)72 and 8772 are tabulated instead of the 
straight differences, the interpolation operation involves 
only two additions and two multiplications. This table con- 
tains ten entries, requires second-order interpolation, and 
the maximum error is 4 in the fourth decimal. A punched 
card table may be key punched directly from a printed table, 
it may be constructed with the use of machines, or it may 
be copied from someone else's punched card table. One 
value of the argument, its function value and interpolating 
aids are punched in a single card so that each card car- 
ries a line of the table. For example, a card for Table III 
would carry x^, sin Xi, (8f_i/2+8i+i/2)72and 87/2. A dis- 
tinguishing X is punched in some column for identification 
purposes. 

Assume that we are given a set of detail cards, each of 
which carries a value of the argument to be used in consult- 
ing a table. The arguments in the detail cards are presum- 
ably the result of a previous calculation. The first step in 
consulting a punched card table is to merge the detail cards 
with the table cards in such a manner that each table card 
is followed by all the detail cards whose arguments are 
greater than or equal to the argument in the table card but 
less than the argument in the following table card. This 
may be done on the sorter (if the arguments in the detail 
and table cards are in the same columns) or on the collator. 
If the sorter is used, the two sets of cards are arranged in 
ascending order of the argument, with the table cards enter- 
ing the machine first. If the collator is available, and each 



set of cards is in ascending order of argument, then a simple 
run will accomplish the same purpose. The table cards are 
placed in the primary feed of the collator, and the detail 
cards are placed in the secondary feed. The argument of 
the first table card is compared with that of the first detail 
card for a high, low, or equal condition. If the arguments 
are equal, both cards are ejected simultaneously, with the 
card from the primary feed alighting first. If the arguments 
are not equal, the card with the lower argument is ejected. 
The control panel is wired to perform these operations. 

In the case of the critical table it is only necessary to 
transfer the value of the function from a table card into the 
following detail cards with the use of the gang punch. The 
distinguishing X punch in the table card prevents punching 
from the last detail card of a group into the following table 
card. If the table is one in which interpolation is necessary, 
the cards may be passed through a calculating punch imme- 
diately after the merging operation. The values of the func- 
tion and the differences for a table card are stored in the 
machine and used for the following detail cards. The next 
table card signals the machine to reset the previous data 
and replace it with new values. 

Some features of the IBM Type 602-A Calculating Punch 
will be discussed, and a description of the control panel 
wiring for second-order interpolation will follow. There are 
30 counter positions divided into six counter groups, and 
six storage registers of 12 positions each. All storage regis- 
ters except the first one are divided into a right- and a 
left-hand group of six digits each. The first storage register 
is divided into an eight-digit right-hand group, which stores 
the multiplier and divisor, and a four-digit left-hand group. 
Two of these storage registers are capable of punching re- 
sults. The counter groups are capable of adding and sub- 
tracting, and do not reset until they receive an impulse to 
do so. The storage units do not accumulate, and they auto- 
matically reset before a new number is accepted. There are 
12 available program steps, which may be repeated or ex- 
tended. In multiplication, the multiplier is read into the 
right-hand part of storage register 1 (IR), and the multi- 
plicand is read out of any storage unit or counter into any 
other counter to accumulate the product. 

The programming for the second-order interpolation, 
U = fi-\-u [(8^^1/2+81+1/2)72 + u (87/2)], is as follows 
(Figure 1): 
Table cards: Read ji (function) into storage unit 3R. 

Read (8i'_i/2+8^+i/2)72 (modified first dif- 
ference) into storage unit 2L. 
Read 8"72 (modified second difference) 
into storage unit 2R. 
Detail cards: Read u into storage unit IR. 



CALCULATING PUNCH - TYPE 602A - CONTROL PANEL 






8 9 10 n 12 13 14 
-CONTROL READING 



I<5 17 18 19 20 21 22 
-PUNCH 




23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 



—7 — 6-PUNCH EMITTER-1 — 11 — 12-1 

ooooooooool 



Figure; 1 
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Program step i: Multiply u 8i'/2: Read out Bi' /2 from 
storage unit 2R, and develop the negative product in count- 
ers 5, 6. 

Program step 2: Add (8Ui/2+8i+i/2)/2toM(8i')/2:Read 
out ( 8i_ 1/2 — 8i+ 1/2)72 from storage unit 2L, and read posi- 
tively into counters 5, 6. 

Program step 3: Multiply [ (8i_i/2+8i+i/2)/2 + u (8'/)/2] 
by u: Read out [(8i'_i/2-f-8Ui/2)/2 -f u {Si')/2] from 
counters 5, 6 and develop the product in counters 1, 2, 3. 

Program step 4: Add/i to u [ (8^1/2-^8^^+1/2) /2 -|- u (81') /2] 
to obtain /„: Read out fi from storage unit 3R and read posi- 
tively into counters 1, 2, 3. Reset counters 5, 6. 

Program step 5; Punch the result /„• Read /„ out of counters 
1, 2, 3 into punching storage unit 6. Reset counters 1,2, 3, 
reset to 5. 

The design of a punched card table for a particular prob- 
lem depends upon many factors. If the problem is of con- 
siderable magnitude, it is worthwhile to spend time and 
energy to prepare special tables. On the other hand, for a 
shorter problem a less refined table may be more efficient. 
In general, the selection and design of the punched card 
table depends upon the problem, the machine operators, and 
the machines available. Each table should be constructed so 
that it will be most economical. 

Usually, it is necessary to trade table size for length of 
computation. Consider, for example, Newton's formula 



/. 



f, + uM+'<^Ai'+ 



If a critical table is desired, the effect of the first differences 
must be negligible; therefore, all possible function values 
are needed. If only first differences are permitted, the effect 
of second differences must be negligible; however, a larger 
interval than that in the critical table is permitted, and 
fewer entries will be necessary in the table. The computa- 
tion of the function now involves an addition and a multi- 
plication. In the three tables we discussed previously, the 
critical table contained 1001 entries, and its accuracy was 
0.5 X 10~^; the table with linear interpolation contained 91 
entries and had the same accuracy; the table with second- 
order interpolation contained 10 entries, and its accuracy 
was 0.4 X 10""'''. This process of the reduction of table size 
can be continued until one table entry remains, but the 
amount of necessary calculation is huge. In fact, in many 
cases this calculation may be equivalent to evaluating the 
function by its series expansion. In constructing a table, it 
should be kept in mind that the first differences vary 
directly with the interval, the second differences vary di- 
rectly with the square of the interval, etc. 



In most printed tables, the tabular interval is a unit in a 
given position of the argument, and it is usually constant 
throughout the table or varies by a power of 10. This choice 
of interval facilitates the interpolation process for the com- 
puter, and hence fewer errors result. In many problems it 
is desirable to use the largest tabular interval in which it is 
legitimate to interpolate with a given order of differences. 
A table constructed on this principle is called an optimum 
interval table. The punched card method conveniently han- 
dles varying intervals, and thus makes possible a consider- 
able saving in the number of entries in the table. Moreover, 
the formula used for «th order interpolation in the optimum 
interval method is equivalent to that used for a constant 
interval table. 

Consider a linear optimum interval table of the exponen- 
tial function, y = e'", in the range .s; = to .r = 5, where 
the intervals are so chosen that the error due to neglecting 
second differences does not exceed .05. The permissible 
intervals have been computed, and a portion of the table 
follows: 



X, 



Interval 



P 



D 



1.06 
1.37 
1.64 



31 
27 



.70 



3.38 



4.85 
4.90 
4.95 



It is clear that a large interval may be chosen at the begin- 
ning of the table, but the interval must decrease sharply 
near the end of the table. To facilitate the interpolation 
operation, the formula for ordinary linear interpolation is 
transformed so that, 

f (^i+i) - /(^i) 



e'" = P -\- Dx where D = 



^t+i — ^i 



and P — j{xi) — XiD . 

For example: 

^1-23 = -.70+ (1.23) (3.38) = 3.46. 

The true value for e^-^^ is 3.42, which is within the limit of 
accuracy. This table consists of 40 cards; a table with uni- 
form intervals of .01 contains 500 cards. It should be noted 
that the machine operations are the same as those necessary 
for ordinary linear interpolation. 
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Consider the problem of determining the sines of cer- 
tain angles in the first quadrant with a given accuracy of 
1 X 10"'^. The chart below gives the type of table, the order 
of interpolation, and the number of necessary cards. 





Order of 


Number 


Type of Table 


Interpolation 


of Cards 


critical table 


first order 


5,000,000 


constant interval 


first order 


9,000 


optinjum interval 


first order 


1,500 


optimum interval 


second order 


60 


optimum interval 


third order 


15 



It is apparent that for different problems different tables 
are more efficient. The third-order optimum interval table 
containing 15 cards looks attractive at first glance, but if it 
is to be used only a few times, it would not be economical 
to construct. 

re;ference;s 

1. W. J. EcKERT, Punched Card Methods in Scientific Computa- 
tion. Revised edition in preparation. 

2. H. R. J. Grosch, "The Use of Optimum Interval Mathematical 
Tables," Proceedings^ Scientific Computation Forum, 1948 
(IBM), pp. 23-27. 

DISCUSSION 

Dr. Sherman: In the problem of inverse interpolation, 
can optimum interval tables be used ? 

Miss Krawitz: It is not possible to use the optimum in- 
terval table directly for inverse interpolation. The values of 
F and D must be modified. 

Mr. Opler: There is an article in the July MTAC on 
maximum interval tables. They appeared to be differentiat- 
ing them from optimum interval tables. 

Dr. Grosch: Stadler's article, which was the first sent to 
this country, is just an example of the confusion into which 
we are always getting. He prefers the term maximum inter- 
val, because if you apply the methods which he worked out 
independently of me for determining intervals, the calcula- 
tion of the maximum allowable interval is obtained, but 
when Clemmons and Herget coined the phrase optimum 
interval, they implied a little more than choosing the largest 
possible interval. They implied the possibility of, having 
found the maximum interval, abandoning it in favor of a 
more suitable interval, if something is to be gained by it. 
In the list of tables of the Watson Laboratory there is a 
third-order uniform interval table of sines and cosines 
which is accurate to about 8^ decimal places. That is a 
uniform interval table of 25 lines per quadrant, because, both 
sine and cosine were present in the same line of the table, 
and the allowable maximum interval governed, for half of 



the quadrant, the maximum function allowable for the co- 
sine function. Actually, we would have only reduced the 
table to about 23 lines per quadrant if we had taken the 
maximum interval table. 

We gained all sorts of advantages in implementing the 
table on the mechanical computing aids available. I think 
we ought to say that an optimum interval table is a table 
best adapted to the facilities at your disposal. That is, per- 
haps, the initial stage that you go through in arriving at an 
optimum table. We have many tables at the Watson Labo- 
ratory which started out to be extremely elaborate ones and 
finished in uniform interval form because we could use / 
instead of P and make the table more easily recognizable 
and usable. 

As I see it, all phases of the problem, choosing the proper 
interpolation, forming the proper intervals of the table, and 
so forth, must be considered in order to make the table as 
economical as possible to use for the particular problem you 
are trying to solve. 

Lt. Hastings: I would like to ask for some general com- 
ments on the extent to which it is no longer necessary to 
have tables of simple functions like sine and cosine, because 
of the program repeat on the CPC. 

Dr. Grosch: There are certain mathematical functions 
that are very well adapted to the idea of direct mathematical 
handling. I think an example is the sine, cosine, where you 
have a very simple formula which converges very rapidly, 
for which only every other power of the power series in- 
volved is required, and in which all of the coefficients are 
simple rational numbers. Obviously, you can use the pro- 
gram repeat device and evaluate the sine function very 
nicely without reference to any of the tables discussed. The 
critical idea, in my mind, is whether or not you are able to 
step from value to value. If your problem is — such as it 
might very well be in a trajectory, for instance — that each 
table lookup is the value of the argument which is pretty 
near the preceding one, then you add the important ad- 
vantage of being able to use, what I might think of as your 
past value of the X function, for instance, and obtain the 
X -f- AX, using AX as a variable. In some cases you must 
be very careful about accumulation of errors. On the other 
hand if you have random lookups, evaluation of each func- 
tion, without the aid of a table, is more advantageous. 

Dr. Hurd: I think what prompted Lt. Hastings' question 
is the fact that he has made a survey of computing facilities 
around the country. He noticed in many places where a 
CPC is available that such functions as trigonometric func- 
tions and square root, etc., were being calculated directly 
on the calculator as these functions were needed. 
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SIMULTANEOUS linear algebraic equations up to 
and including order 21 may be readily solved on the IBM 
Card-Programmed Electronic Calculator by using a basic 
approach sometimes referred to as the Crout method.^ A 
slight modification has been incorporated in the procedure 
in completing the back solution. Having obtained the aux- 
iliary matrix at the end of the forward solution^ the operator 
rearranges the elements in such a way that almost the same 
reduction technique may be used for the back solution. 

The method involves the augmentation of the matrix of 
coefficients of the unknowns. To this original matrix is 
added one of the following: . 

a. One column of constants, if a single solution is de- 
sired; or 

b. Several columns of constants, for several solutions; or 

c. A unit matrix, if the inverse is desired. 
Regardless of which end result is needed, the reduction 
calculations are the same. 

By using the CPC, the solution is obtained much more 
rapidly and with considerably less card handling and proc- 
essing than has been possible heretofore on IBM calculating 
punches. 

Renumbering Elements of Original Matrix 

To take advantage of the internal storage numbering and 
to facilitate coding instructions, the familiar notation for 
each element of the matrix has been altered. Instead of 
using 1, 2, 3, .... to designate row i or column j location, 
one numbers the elements after the code numbers of storage 
units and counters in the CPC. 

It takes six hours to invert a 20 X 20 matrix on the 604, 
sorter, reproducer and accounting machine, and it takes two 
hours to perform the same inversion on the CPC. 
The storage and counter designations are as follows: 
Storage bank 1: 11, 12, 13, 14, 15, 16, 17, 18 
Storage bank 2: 21, 22, 23, 24, 25, 26, 27, 28 
Counter groups: — , 2, 3, 4, 5, 6, 7 



The code digit 7 before any counter group number is, 
here, an instruction code meaning add into. Thus, it is ad- 
visable to assign row and column numbers as follows: 

Usual notation: 123456789 10 11 
12 13 1415 1617 18 19 2021 

New notation: 11 12 13 14 15 16 17 18 21 22 23 

24 25 26 27 28 72 73 74 75 76 

Counter group 7 is reserved for the accumulation of the 
check sum. Counter group 1 is not used to store element 
values. Thus, it is possible to store 21 elements and provide 
for the accumulation of a check column. Each counter- 
group and storage location has the capacity for handling a 
ten-decimal digit number and its associated sign. 

Augmenting the Original Matrix 

Before the reduction calculations are made, the original 
matrix, consisting of the coefficients of the unknowns, must 
be augmented. The augmentation is always made by adding 
columns of elements to the right of the original matrix. 

To facilitate identification, the composite matrix, includ- 
ing the original matrix, the augmented portion, and the 
check column, is divided into sections called subdivisions. 
The original matrix is subdivision 1, the augmented portion 
is subdivision 2, and the check column, each element of 
which consists of the sum of the elements in the correspond- 
ing row, is subdivision 3. The augmented portion, subdivi- 
sion 2, will vary, depending on whether one or several solu- 
tions are desired, or whether the inverse matrix is to be 
obtained. In this paper the case of the inverse will be treated. 
The other cases would be handled as described below, but 
the augmented portion would consist of as many columns 
of constants as the number of solutions required — instead 
of the unit matrix. 

Using the new notation and augmenting the original 
matrix so as to prepare for calculating the inverse, one has 
the twenty-first order matrix shown in Figure 1, page 58. 
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Originai, Matrix 
Subdivision: 1 

^11,11 ^11,12 0^11,13 . . . ail, 75 %1,76 
^12,11 ^12,12 ^12,13 • • • ^12,75 ^12,76 



%5,11 ^75,12 ^75, 13 . 
^76, 11 ^76,12 ^76, 13 • 



<^75,75 ^76,76 
<*76,75 ^76. 76 



ill, 11 
0l2,ll 



075,11 
076,11 



Unit Matrix 

2 

Oil, 12 • • • Oil, 75 

ll2,12 • • • 0i2,75 



^^75, 12 
076,12 



0ll,76 
0i2,76 



175,75 075,76 
076,75 176,76 



Che;ck 

Column 

3 

Cii,ii 

Cl2,ll 



<^75,11 
^76,11 



Figure 1 



1. Take the corresponding element aij in the composite 
matrix. 

2. From it subtract the products resulting from pairing 
the elements of the auxiliary matrix in the row to the left 
and in the column above the desired element, taking the 
products in order, i.e., first in row by first in column, sec- 
ond in row by second in column, etc. (This rule will be 
illustrated later.) 

3. If the desired element is on or below the principal 
diagonal, record the result obtained in 2. If the desired 
element is to the right of the diagonal, divide the result 
obtained in 2 by the diagonal element in the same row of 
the auxiliary matrix. 

The reduction calculations are applied to each element of 
one row of the original composite matrix to obtain the 
corresponding row of the auxiliary matrix, completing the 
solution for that row for the forward solution. Calculations 
begin on the first row and continue on succeeding rows 
until each row of the original matrix has been reduced to 
form the auxiliary matrix. 

After one has developed a portion of the auxiliary matrix, 
it will appear as in Figure 2. 

The zero elements to the right of the diagonal in the unit 
matrix need not be introduced until the back solution is 
started. 

Reduction Calculations 

To obtain the inverse matrix (or specific solutions) each 
element of the original composite matrix is reduced by a 



series of calculations which may be divided into two major 
portions. The first portion consists of the reduction calcu- 
lations applied to each element of the original composite 
matrix. These reduce it to a matrix known as the auxiliary 
matrix. This first reduction is called the forward solution. 
The second portion of computations applied to each element 
of the auxiliary matrix (whose elements have been shifted 
as described later) reduce it to the final result which is 
the inverse matrix. This second reduction corresponds to 
the back solution of the Crout method. 

Regardless of whether the reduction calculations are for 
the forward solution or the back solution they follow the 
same basic rules. These rules for calculating the value of a 
reduced element bij are: 

As an example of how to apply the rules for calculating 
the value of a reduced element bij of the auxiliary matrix, 
consider the computation 

^14,15 = 

%4,15 + (~^14,ll ^11,15 ~ ^14,12 ^12,15 ~ ^14,13 ^13,15) 
^14,14 

During the calculation of the reduced elements bij in any 
one row of the auxiliary matrix, the elements to the left 
and on the principal diagonal are stored, and all of the 
elements to the right of the principal diagonal are summary 
punched as they are calculated. Each of these summary 
punched cards becomes the source of one of the pair of 
factors involved in negative multiplication. At the comple- 
tion of the calculation for one row, the summary punched 
cards for that row are merged in with cards summary 



Subdivision: 1 

'*-i'll,ll'^ -^41^12 ^11,13 ^11,14 &11,15 • • • 

'^12,11 ''■>^12,12 "■-^2,13 ^12,14 &12,15 • • • 

^13,11 ^13,12" ^ i'lS.lS ^ -^13,14 ^13,15 • • • 
^14,11 ^14,12 &14,13->^14,14"- 

"HPrincipal diagonal 



^11,76 


^11,11 




^12,76 


^12,11 


^12,12 


^13,76 


1 ^13,11 


^13,12 



'13,13 



3 
^11,11 
C?12,ll 
C^13,ll 



Figure 2 
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Subdivision: 

^11,12 &11.13 
^12,13 



1 

^11,14 
^12,14 



• ^11.75 

• ^12,76 



^74,75 



^11,76 
^12.76 



^74,76 
^76,75 



^11,11 
&12,11 



^r4,ll ^74,12 
^7»,11 ^76,12 
^76,11 ^76,12 



I 
I 

Figure; 3 



^74,74 

^76,74 ^76,75 

^76,74 &7e,76 



'76,76 



3 

^11.11 
C^12,ll 



^74,11 
0^76,11 
^76,11 



punched from previous rows. These combined cards are 
arranged in groups to form the columns of elements above 
the particular row elements to be reduced by the next cal- 
culation. The original elements ai,^- are also in this deck. 

As the calculations proceed across a row, each result to 
the left of and on the principal diagonal is stored internally 
in the storage location corresponding to the identification of 
the column of the original element involved. In the example 
above, the &i4,i5 element will be stored in the column desig- 
nation of ai4,i5 which is storage 15. Similarly, code76means 
add into counter group 6. 

While the summary punched cards in the column above 
the particular element being calculated feed through the 
IBM Type 402 (or 417) Accounting Machine, each card is 
paired with its proper other factor for negative multiplica- 
tion. This proper pairing is accomplished by having that 
storage location designated by the column number read-out 
to combine with the summary punched card containing the 
same row number. That is, the summary punched card con- 
taining one factor hu,j and row designation k calls out the 
second factor hi^it from the internal storage of the same 
column number k. 

As the original element card feeds through the type 402 
accounting machine, its value is added. 

For any one row the divisor is always the same and is 
available once the principal diagonal element has been cal- 
culated. The storage location of the divisor is remembered 
for the entire row and is used to complete division at the 
proper time. 

The same reduction operation applies throughout the en- 
tire composite matrix, including the unit matrix and check 



column. At the end of the forward solution the auxiliary 
matrix will be as shown in Figure 3. 

To identify correctly, by row and column, results that 
will appear in the spaces occupied by the zero elements in 
subdivision 2, the following cards must be added: 

Oil, 12 Oil, 13 • • • Oil, 75 Oil, 76 
0l2,13 • • • 0i2,75 0i2,76 



074,75 074,76 
075,76 

Back Solution 

The Crout method back solution is modified so that the 
same calculating rules may be used as were employed in the 
forward solution, except that no division is involved at any 
time. To do this it is necessary to shift the arrangement of 
the elements in the auxiliary matrix before starting the re- 
duction calculations. 

This necessary rearrangement of elements is rapidly car- 
ried out on the sorting machine. The necessary shifts are: 

1. Invert the order of the rows in each column of each 
subdivision, at the same time combining the zero elements 
with subdivision 2. 

2. Invert the order of the columns in each row of sub- 
division 1 and 2. 

Thus, the rearranged auxiliary matrix is as shown in 
Figure 4. 



Subdivision: 

^76,76 

^74,76 ^74,75 



^12, 75 
^11,75 



1 



l^l2,14^^ 
^11.14 



bl2,13 
^11,13 



j ^76, 76 
.O7 



^11. 



''75,73 







11,76 



^76,75 
^75.75 







11,75 



Figure 4 
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^76,12 


^76, 11 


<^76,11 


^75,13 


^75,12 


^75,11 


C^75,ll 


^74,13 


^74,12 


^74,11 


<^74,11 




&12,12 


^12.11 


<^12,11 




Oil, 12 


&11,11 


rfll,ll 
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Repeating the calculations except for division, one obtains 
the rearranged auxiliary matrix reduced to its final form. 
Subdivisions 2 and 3 are of interest here: 



Subdivision 

I 



I ^76,76 ^76,75 
I ^75,76 ^75,75 



I ^12.78 ^12,75 

I 

I ^11,76 ^11,75 



^76.12 ^78,11 j /76.II 

^75,12 ^75,11 I /76,11 

I 



I r 

^12,12 ^12,11 I /12,11 

I , 

^11,12 ^11,11 I Jll.H 



Subdivision 2 is now the inverse matrix. The inverse 
matrix may, of course, be used in matrix by vector multi- 
plication with any column of constants to obtain the specific 
solution corresponding to those particular constants. If one, 
or several, columns of constants had been used in the aug- 
mented portion of the matrix, the specific solutions would 
appear in as many columns of subdivision 2 above in place 
of the inverse matrix. 

The idea of using the Crout method on the CPC was 
originated by Mr. William W. Woodbury, Northrop Avia- 
tion Company. 

Details of planning charts and wiring diagrams may be 
obtained from the Applied Science Department, Interna- 
tional Business Machines Corporation, 590 Madison Ave- 
nue, New York 22, New York. 

1. WiLUAM Edmund MunE, Numerical Calculus (Princeton Uni- 
versity Press, 1949), pp. 17-25. 

DISCUSSION 

Mr. Blkins: Is this method the most efficient method now 
known for inverting a matrix ? 

Dr. Hurd: This is the most efficient method that we have 
so far discovered for the card-programmed calculator. As 
far as we know, this minimizes the amount of summary 
punching which is necessary. On standard machines, using 
the combination of any one of the calculating punches, the 
sorter and the reproducer particularly, if there are a num- 
ber of matrices of the same order to invert, you will also 
have very efficient operation. 

Dr. Sherman: There are several things I might say. If 
the original matrix is s)mimetrical, which is true for the 
case of normal equations, then the auxiliary matrix can be 
calculated slightly more easily because the elements below 
the diagonal are related to the elements to the left and right. 
You can either save almost half of the calculations, if you 
have a case of a symmetric matrix, or you can carry out 
this procedure and have an additional check. There is one 
obvious difficulty which actually arises; that is, if any one 



diagonal element is very small or close to zero, there are 
difficulties with the later calculations. How do you handle 
that case? 

Dr. Grosch: That is a hard problem. I inverted a 19 by 19 
where I knew in advance I was going to have a lot 
of indeterminates. I used the 602-A punch and carried 
24 X 24 = 24 multiplications throughout. I lost 16 figures, 
which were so many that I could not iterate by the simple 
iteration procedure. I must do it over again with 32 figures, 
one of these days. There are other methods of using what I 
obtained. The trouble is that we have all sorts of interesting 
ideas about how to handle the problem of troublesome mat- 
rices but usually we don't know they are troublesome until 
we have started the problem. My own idea is that the 
Goldstein estimates of the laws of accuracy in inverting 
large matrices need not bother us unless the matrices arise 
from something like operational research where you have 
no idea at all of how the matrix is going to behave. 

The kind of matrices that arise in flutter work are not 
likely to give trouble unless, in the few cases such as those 
with which I was experimenting, you have been warned in 
advance by the physical situation. I was using this as a 
matrix for fitting polynomials to curves which did not fit 
well, and I knew that determination of the coefficients was 
going to be hard; so I was not surprised at my troubles. 

All of our discussion has been matrix arithmetic so far 
and has centered on the triangulation theme. There will be 
another paper on the Seidel or relaxation procedure. There 
is one other process worth considering on the card-pro- 
grammed calculator. If you have a symmetrical matrix with 
which to deal and want an inverse in a simple manner, you 
may use the square root method which is now known under 
the name of Choleski. This is a process of finding a tri- 
angular matrix which, when multiplied by its transpose, 
equals the original given matrix. Having a triangular ma- 
trix you can find its inverse by operations similar to the 
back solution in the Gaussian elimination scheme. 

This may be well worth considering, but the job of find- 
ing the triangular matrix is hard on any kind of standard 
equipment. With the card-programmed calculator, a pack 
of eight or ten cards is needed for a second-order matrix, 
then with a few more added for a third-order matrix and 
fourth-order, etc., up to several drawers full for 20th order. 
You would invert your matrix by loading the whole thing 
into the card-programmed calculator, and summary punch 
the results without any wasted time summary punching. 
Then, the cards just summary punched would be placed in 
the machine for the back solution, which is very simple. 

Mr. Opler: I would like to make a general remark also 
about matrix inversion. Because the usual ehmination pro- 
cedures have so many different steps, I have spent consid- 
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arable time investigating methods of matrix inversion other 
than the Gaussian ehmination or row-by-row method. 
There is no doubt that this method which you have demon- 
strated is the only method that will work for virtually all 
cases of matrices. However, at times you have a borderline 
case which may be troublesome if yoti are not careful. There 
are a number of methods that appeal to me because I have 
done my calculating through the accounting department per- 
sonnel, and I find that they have difficulty in a long proce- 
dure involving row-by-row elimination; so I prefer, in gen- 
eral, a method that will work by a simple procedure. There 
are several simple procedures for inverting matrices, but 
they do not apply to all matrices. One, which is without a 
doubt the simplest method and which I have carried out, is 



the Monte Carlo Method. The entire method consists in 
taking a pack of cards and putting them in an accounting 
machine. It's a very nice method, but unfortunately is only 
adapted to classes of matrices whose largest characteristic 
roots are less than 1. 

Mr. Chancellor: I would like to second the motion of 
getting your operation as completely automatic as possible. 
That is the whole IBM goal in all of your applications — 
mathematical, accounting, or otherwise. 

Lt. Hastings: I just wanted to remark that using the 
card-programmed calculator simultaneously for two sys- 
tems, because of the alternation of the calculation and the 
collating steps, two systems could probably be done in very 
nearly the same time as one system. 
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The Gauss-Seidel NLethod of Solution 
of Simultaneous Linear Equations 

THE GAUSS-SEIDEL method for the solution of 
simultaneous linear equations was first advanced by Gauss 
and Seidel about the turn of the 19th century in connection 
with the solution of normal equations. In addition to the 
solving of normal equations, this method may be used to 
solve various other systems if they meet certain require- 
ments. We have used this method successfully, and the air- 
craft industry uses it very extensively, finding that it yields 
good results. 

I would like to list the advantages of this method in an 
attempt to show you, as we go through the procedure, how 
these things are advantageous. We can obtain accurate an- 
swers to a large number of digits at a rapid rate. The prob- 
lem is easy to set up for the machine. The machine runs 
continuously until an answer is obtained. In one 10 by 10 
set of equations considered, results were obtained in thirty 
minutes, with the answers to seven decimal places. In about 
twenty minutes, the answers to five places were obtained. 
In less time than that, there were three or four places. The 
operator is in touch with the problem constantly and can 
tell how the problem is progressing immediately, since each 
new value of the variable and the correction factors are 
listed after they are calculated. The operations are com- 
pletely repetitive, and almost any standard control panel 
can be used. 

The mathematics is quite simple, and I would like to go 
through it briefly with you. Consider a set of three algebraic 
equations, bearing in mind that this method is general for 
any number of equations: 

axxXx + ai2^2 -|- ais-^'a - &i = , 

«21-^l -f ^22^2 -|- CLlzXz — &2 = , 
a^lXx -f az2^2 -\- «33^3 — &3 = , 

Each of the elements, a^,;, is on a separate card. If the values 
of the x's are known, each of these equations will yield zero. 



to some accuracy. Assuming these values are unknown, let 
a first guess be made for each unknown x. Now the equa- 
tions will not yield the answer zero, but something different 
from zero, say e^, and the equations will be as follows: 

ailA'i -|- ai2^2 + 013-^3 - &i = ei , (1) 

ff21-«'l -f- a22X2 + ^23^3 " &2 = €2 , (2) 

«31-^1 + a32-«'2 -f- a 33-^3 " &3 = ^3 • (3) 

The general procedure is this: 

1. Using the first guess for Xx, X2, and Xz, solve for ei. 

2. Determine a new Xx that will cause this equation to 
yield zero from this formula 

€ 

^1 new ^^ ^1 old "7 > 

da 

where e = residual from an equation and an is a main 
diagonal coefficient of that equation. For equation ( 1 ) 
a is axx and for (2) a is a22, etc. 

3. The new Xx is used in place of the old Xx from now 
on until a better value is found on the next cycle. 

4. Using the new Xx, the old X2 and xz, solve for €2 from 
equation (2). 

5. Solve for the new X2 to make the equation go to zero. 



.STrl 



•^old ~" 



£2 

0^22 



6. 
7. 



9. 



Substitute the new X2 for the old x^. 

Solve for es from equation (3). 

Solve for the new x^. 

Substitute the new x^ for the old Xz. This is the end 

of a cycle for a set of three equations. 
The process is repeated using one new value each time 
an equation is solved. As this process is continued, the val- 
ues will converge to the correct answers for all equations, 
and the residual values will approach zero. This method has 
been found to converge for nearly all algebraic equations 
derived from engineering and physical problems and normal 
equations. 
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Machine; Procedure; 

Assume that a standard 604 calculator control panel is 
available that will do at least the following operations nec- 
essary for this problem: multiply, divide, subtract, and 
transfer from channel C to channel B. 

The 402 accounting machine control panel is the standard 
control panel that reads factors from channels A and B and 
coding instructions from the card. It is wired to iviST from 
from an X in the OP field channel A,B,C and instructions. 

The additional wiring consists of the following: all of the 
ttij- coefficients are in cards; as the main diagonal, au, is used 
to multiply by the proper x, it also must be stored to be used 
in the calculations for the new x. 

€ 
Clii 

Channel A is wired to the entry hubs of counter group 7. 
An identifying X-40 is punched in all cards containing main 
diagonal, an, coefficients. This X is used to pick up pilot 
selector 1 1 from the upper brushes. The sign of coefficient a 
will pick up pilot selector 12. When the main diagonal co- 
efficient is read from the lower brushes, it will be added or 
subtracted into counter group 7 by the following wiring: 



No. No. 

11 12 

T ^ • — »-Group 7 Subtract 

Channel A 



oooooooooo 
o — e — e — e — g — e — e — o — e — o 




Card Cycles Pulse 



Counter Entry Group 7 

(could be selected with two coselectors) 



FiGURK 1 

Channel A is also wired to the entry hubs of counter 
group 7.* 

Coding of Cards 

There are 22 ten-position storage registers in the CPC. 
Sixteen of these are registers, and six are counter groups. 
Let one counter group, say counter group 6, be used to ac- 
cumulate the sum of products, ^ax — b = e. Group 7 is 
used to store the main diagonal coefficient an. That leaves 
20 storage units for as many as 20 unknown jit's. Each x is 
assigned to a storage unit and will be called for use from 
that unit at the proper time. When a new x is developed it 

»With this added wiring, counter group 7 cannot be read out with- 
out resetting on channel A unless the wiring from channel A to 
the entry hubs is put through the transferred points of the co- 
selectors coupled to pilot selector 11. This is unnecessary in this 
application as counter group 7 is always read out and reset. 
Counter group 7 cannot be read out on channel B at all; it will 
not read out properly. 



will be sent to the storage unit to replace the proper old x 
value. Xi and X2 will be assigned to counter groups 2 and 3, 
respectively. All other values may be assigned in either' 
counters or storage registers. This is done to avoid un- 
necessary blank cards. All of the coefficients, the values of 
the ay's and the right-hand column, or b/s, are key punched. 
There will be as many cards as there are Oi/s and b/s, one 
value on a card. Each card will also contain row and column 
identification of the coefficient as it appears in the matrix 
form, and instructions to the machine to perform the re- 
quired operations. 

Consider equation ( 1 ) of a 20 X 20 set of equations: 

^ll-*"! + ai2X2 -\- ais^'3 + aiiX4 -f- . . . -h «120-*"20 — bi — ei . 

Card I contains an which is read on channel A. It con- 
tains instructions to read Xi on channel B and multiply Xi 
by an, adding the product in counter group 6. 

Cards 2 through 20 contain similar instructions. Each 
card contains an c^y coefficient and calls out the proper x, 
multiplies, and accumulates the answer in counter group 6. 

Card 21 reads the right-hand side through the calculator 
to subtract in counter group 6. 

Card 22 is a blank card necessary because the accumu- 
lated answer e is wanted for the next operation, and it is 
not complete until the right-hand side &_,• is subtracted. 

Card 22 reads out e which is in counter group 6 on chan- 
nel A. It reads out the main diagonal coefficient, a«, which 
is in counter group 7, and divides c/aii, and transfers the 
answer back to channel B. This answer will be listed from 
channel C. 

Card 2/\. reads out the proper x^v^ on channel A and per- 
forms the operation 



•^old 



— Xt, 



.^■new is stored in place of x^xa ready for use from then on. 
There are a total of 481 cards. 

The cards for the next row initiate the same operations 
of multiplying and accumulating the answer in counter 
group 6. The main diagonal, a^i, is stored automatically and 
used in card 23 of each row to calculate (./an. Card 24 cal- 
culates the new value of the x concerned and stores it in 
place of the old value to be used in the following operations. 

This deck is fed through the machine as many times as 
necessary for desired convergence. For convenience, an- 
other deck may be 80-80 reproduced from the correct deck, 
thus giving 962 cards. Cards may be taken from the card 
stacker and replaced in the feed hopper, thereby keeping 
the machine in constant operation. 

It would be very helpful if a good first guess could be 
entered in the machine. A good first guess would be the 
first digit of the actual answer. The machine would imme- 
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diately show refinement of these answers. This in general 
is not necessary, fortunately, and the method usually works 
nicely with an initial guess of zero for all values. This is 
done by resetting all the storage registers and counters to 
zero in the beginning of the process. 

One can tell that the process is functioning properly by 
watching the printed indications of the correction factors 
and the values of the x. As the process converges, the cor- 
rection will become smaller and the x values will remain 
constant. In this process the machine may be run non-net 
balance to avoid conversion cycles. All answers are listed 
from channel C and are true figures with a sign. 

If a master deck is available for the 20 X 20 set of equa- 
tions, it will take approximately 2^ hours to key punch the 
aij and hj values into the cards, prepare the deck by repro- 
ducing, and run the problem on the CPC. An example of 
the solution of a typical 20 X 20 after it was started on the 
CPC is as follows: A first guess of all values oi x = is 
made, with eight-decimal digits in original data and seven- 
decimal digits of accuracy required in the results. Twenty- 
seven iterations (times through the entire set of twenty 
equations) are completed in 90 minutes. Complete setup 
time, including punching of original data, should be less 
than one hour, using efficient methods. 

Recommended Setup Method 

A master deck for the solutions of a 20 X 20 set of equa- 
tions is on file. It contains all necessary punching, except 
that the field in the cards set aside for the coefficients to be 
read on channel A will be blank. Assume that a set of 
14 X 14 equations is to be solved: 

The master deck is in order for operation on the CPC, 
that is, row, column. There are only 14 equations, so that 
the cards for row 15-20 will be taken out by hand and set 
aside. Sort the master deck in column order. This will put 
the card numbers for each row together. Recall that there 
are 24 cards for the solution of one equation in a set of 

20 X 20 equations. Therefore, in column order, or card 
number in the row, there are 24 groups. Pull groups 15-20, 
as there is no use for these for a s6t of 14 X 14 equations. 
The remaining cards are the necessary cards for solution. 
Reproduce this deck (80-80). This deck is called the work 
deck. Put the master deck in order and file. 

Key punch the aij coefficients and the right-hand sides, bj, 
each in a card with column and row identification. There 
will be 210 cards punched in approximately 14 columns of 
the card. This punching must be correct; therefore, verify 
the punching ! Sort the key punched cards in the same order 
as the work deck, column, and row. Take groups 1-14 and 

21 out of the work deck. Reproduce the key punched co- 
efficients into this part of the work deck, comparing the row 



and column for correct card punching. Place this deck back 
with the remainder of the work deck and sort into row and 
column order. This deck is ready for the solution on the 
CPC. It may be reproduced (80-80) several times to obtain 
a convenient number of cards if desired. 

CPC Operation^ 

Reset all storage to zero. Commence running the cards 
with the setup switch on list, each card for one complete 
cycle. Check all 210 punched numbers against the original 
data. The machine will list from channel A. Switch off the 
accounting machine from this point, listing only card 23 
and 24 of each row. 

DISCUSSION 

Mr. Liggett: If the matrix is positive definite, the Gauss- 
Seidel method may be used for the solution of the simul- 
taneous linear equations. The process for discovering 
whether the matrix is positive definite is so laborious that 
it is not worth while. However, if the values on the main 
diagonal are large in comparison with the remainder of the 
elements of the matrix, the method is worth attempting. 

Dr. Sherman: I have used the Gauss- Seidel method and 
found that it is very satisfactory for my purposes. In one 
case, in connection with infrared spectroscopy, I wished to 
solve a set of simultaneous equations with quadratic as well 
as linear terms. The Gauss-Seidel method may be applied 
to this problem. I neglected the quadratic terms and used 
only the linear terms for the first approximation. Having 
obtained the first approximation, I used it to solve the 
quadratic term. By this method the iteration may be con- 
tinued and satisfactory results obtained. 

Mr. Liggett: The problem I have discussed with you was 
a 10 X 10 matrix. It is part of a 20 X 20 matrix which was 
actually solved in our New York Technical Computing Bu- 
reau. The elements of the main diagonal consist of values 
between .6 and .9. The off-diagonal elements in some cases 
were as high as .7, with the major portion of the values 
smaller. Our problem of the 20 X 20 matrix was solved 
with 7 decimal-digit accuracy in 28 iterations. The 10 X 10 
portion which I have used required 28 or 29 iterations to 
obtain the same accuracy. This indicates that the number 
of unknowns does not necessarily affect the rate of con- 
vergence of the iterative procedure. 

bThis method may be extended conveniently to 21 unknowns if the 
sums of products are accumulated electronically. Additional stor- 
age units allow the increase of the number of unknowns. If the 
sign of the answer is known, more than one element may be stored 
in the same unit, thus doubling or tripling the number of un- 
knowns that can be accommodated. 
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One problem which has been solved on our SSEC con- 
tained 168 equations and required only nine iterations. Be- 
cause this method does not take zero elements into consid- 
eration, these cards may be removed from the deck before 
calculation. 

The Gauss-Seidel method may be adapted to the IBM 
Type 604 and 602-A Calculating Punches. Because there is 
not the added storage in these machines, it is necessary to 
punch intermediate results on a trailer card. These inter- 
mediate results are reproduced into another field of the 
card, the wires changed, and the iteration procedure re- 
peated. 

If a minor mistake is made during the calculation, there 
will not be any effect on the convergence, but if a major 
mistake is made, the effect would be the same as though the 
process had been begun again. 

If, after a number of iterations, the values for the Xis 
become increasingly large, convergence will not be accom- 
plished. In general, after each iteration the values of the .^'s 
will improve, but not necessarily every value, every time. 

Dr. Sherman: The Gauss-Seidel method has a more gen- 
eral application than the solution of simultaneous linear 
equations. It may be used for many types of simultaneous 
nonlinear equations. 

Mr. Opler: The Fisher method can be used as an approxi- 
mate test as to whether the matrix converges, that is, if 
cbijCLji < diiCiji, the Gauss-Seidel method will converge. This 
is a sufficient condition, but not a necessary one. 
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Approximating the Roots 
of a Polynomial Equation 

The probi^em of approximating the roots of a polynomial 
equation arises frequently in industrial research. For ex- 
ample, in the aircraft industry it is desirable to calculate 
the roots of an equation which is obtained by equating the 
flutter determinant to zero. 

A common method of solution of polynomial equations 
employs the Newton method of successive approximation. 
Thus, to approximate a root oi j{x) — 0, one uses the 
iterative procedure described by 



"« + ! 






John Lowe has described a process for applying this proce- 
dure to polynomial equations with complex coefficients 
through the use of the IBM Type 602-A Calculating 
Punch.^ Mr. Lowe uses synthetic division to calculate f(xn) 
and f{xn) in one run through the 602-A, and calculates 
Xn+i in a second run. 

The purpose of this paper is to indicate that the card- 
programmed electronic calculator has also been applied 
successfully to problems of this sort. We have wired a 
general-purpose control panel for complex arithmetic with 
ten-decimal digit accuracy. A total of sixty cards for data 
and programming will load the machine and carry out an 
iteration for a cubic polynomial equation with complex 
coefficients. 

Extensions of the method of Newton to non-polynomial 
equations through the use of the CPC are clearly possible. 
For example, iterative solutions of equations containing 
trigonometric or exponential coefficients are possible if 
appropriate general purpose control panels are used. 
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Matrix by Vector NLultiplication 
on the IBM. Type 602-A Calculating Punch 
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PUNCHED CARD equipment handles matrix arith- 
metic in a very efficient manner. Before giving a machine 
application, I should like to discuss briefly the basic prin- 
ciples of matrix arithmetic. 

A matrix is a rectangular array of numbers; each number 
is referred to as an element of the matrix. Two indices are 
associated with each element: the first denotes the row, and 
the second denotes the column in which the element is 
found. For example, the element aij is located in the ith row 
and yth column of the matrix A. The following matrix A is 
referred to as an w • m matrix, where n (the number of 
rows) may or may not be equal to m (the number of 
columns). 



021 



ai2 
O22 



\0'nX 0'n2 



Cl'im\ 
^2m 



V 



There is a matrix arithmetic that includes addition, sub- 
traction, multiplication, and, in a sense, division. The sum 
of two matrices, A -\- B = C,is obtained by the addition of 
corresponding elements, i.e., ay -f- &y = %. Subtraction is 
defined in a similar manner. 

Consider a matrix A that has n rows and / columns, and 
a matrix B that has k rows and m columns. The product, 
A-B,is defined only if / = ^; that is, the number of columns 
in the matrix A must equal the numbers of rows in ma- 
trix B. 



/an . . . au\ /b^ 

O21 ... (I21 ^21 



\flnl 



ani/ 






\h^x 



I 



\CnX 



Cxrn\ 
C2 



•/ 



The product matrix, C, is an ww matrix. The elements of 
the C matrix are defined in the following manner: 



C'i- 



I 



For example, C12 = aii&i2 + ^12^22 + • • • + diihz - 

If we have a set of twelve equations in twelve unknowns 

«1,1 ^1 + «1,2 ^2 + . • • + «1,12 -^12 = bi 
^2,1 -^1 + Cl'2,2 ^2 + . • • + fl^2,12 '*"12 = &2 



(1) 
fll2,l-*"l 4" ^12,2-*'2 "}"••• + ^12,12-^12 = ^12 j 

it may be represented by a matrix and two single-column 
matrices. (A matrix that consists of a single row or column 
is referred to as a vector.) 

/%,! ^1,2 . . • «i,i2 \ /-n \ /bi \ 

= . , or Ax = B. 

\ai2.1 Cli2,2 • . . «12,12/ \-«"l2/ V12/ (2) 

If the multiplication of the matrices is performed, the result 
will be the set of equations ( 1 ) . To solve for the values of x, 
we must multiply both sides of the matrix equation (2) by 
1/A or A~^. The inverse of a matrix A, or A-^, is defined 
as the matrix that, when multiplied by A, yields the unit 
matrix /. 



AA- 



A-^A = I 



/I 
1 



o\ 



\0 ... 1/ 

A discussion of matrix inversion will be given in a later 
paper. 

Multiplying both sides of equation (2) we have 
A-^Ax = A-^B . 
Since any matrix multiplied by the unit matrix yields the 
matrix itself, 
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or /;iri \ M^ ... 0142 \ /&i \ 



W/ 



or 



\-^12/ \^l'2,l • • • fll2,12/ 

^i = flu&i + ai'2&2 + . . . + ai'i2&i2 , 

where i = 1, , . . , 12, 

We have now developed a problem that involves the 
multiplication of a vector by a matrix. There are many 
possible methods of performing this process on the calculat- 
ing punches; the most efficient method depends upon the 
machines available, the order of the matrix, the number of 
digits in the elements, and the number of necessary niulti- 
pH cations. The following method was chosen for its sim- 
plicity. Assume that we have a vector with 12 elements, and 
a matrix of the 12th order; 

(^1,1 . . . ai,i2) • /5i,i . . . &i,i2 \ 



\*. 



^12,12/ 



The 12 columns in the B matrix are divided into four 
groups of three columns each, in the following manner: 

^1,1 h.2 h.3 



^12,1 &12,2 ^12.3 
blA &1,5 &1,6 
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Each row of three elements is punched into a card; there 
will be four groups of 12 cards each or a total of 48 cards. 
Three columns of each card carry the identification; one col- 
umn denotes the group, and two columns denote the number 
of the card in the group. A group of 12 cards is punched for 
the A matrix, so that each card carries an element of the 
matrix. Two columns are used for the identification of the 
element. The elements of the A matrix are reproduced into 
each of the four groups of 5-matrix cards. The cards are 
now in the following order: 



1 


1 


ai 


&1.1 


K2 


^1,3 


1 


2 


a2 


^2,1 


h.2 


&2,3 



12 ai2 bi 



^2,3 



4 1 ai 61,10 ^1,11 ^1,32 
4 2 a2 &2,io ^2,11 t'2,12 



4 12 ai2 bi 



^2,11 



bi2j 



A trailer card with a distinguishing X punch is inserted 
behind each group of (detail) cards. The sum of the prod- 
ucts of the vector by three columns of the B matrix will be 
punched into a trailer card. The computation on the 602-A 
is outlined below (Figure 1, page 68). 



C'12,4 ^12,5 ^12,6 
bl,7 ^1,8 ^1,9 



^^12,7 &12,8 ^12,9 
^1,10 ^1,11 ^1,12 



Detail Cards: 

Cards are skipped out without punching. 

Read cycle: Read ai into storage IR. 

Read bij, bij+i, bi^j+2 into storage units 
2R, 3R and 4R, respectively. 

PI Multiply to obtain Oibij, aibijj^i, aibij+2- 

Read b{,y out of storage 2R, read in count- 
ers 1, 2. 

Read ^jj+i out of storage 3R, read in 
counters 3, 4. 

Read bij+2 out of storage 4R, read in 
counters 5, 6. 



^12.10 ^12,11 ^12,12 



P2 



Read. 



*_AHwULAiiiN»j KUiN^-n- lift OU^A-LUIN IKUL PANEL 



C3S 



23 24 25 7b 27 28 29 30 31 32 33 34 35 36 37 38 39 

r— RESET CO SELECTORS 

o T o o o o 



o N o o 



40 41 42 43 AA 



o T o o o o 

4 
O N O O O O 




Figure 1. Vector by Matrix Muetipei cation (X card must be placed before group so that counters reset) 



SEMINAR PROCEEDINGS 



69 



Trailer Cards: 

PI Read out and reset counters 1, 2. 

12 

Read y aibij into storage 6 and punch. 

i=l 

Read out and reset counters 3, 4. 

12 

Read y ai&ij+i into storage 6 and punch. 

i=l 

Read out and reset counters 5, 6. 

12 

Read y aibij+2 into storage 7 and punch. 

i=l 

This procedure for multiplying a vector by a matrix is 
quite general. If the order of the matrices increases, the 
number of necessary cards increases, but the control panel 
wiring and the process remains the same. If the number of 
digits in each element is too large for the sum of the prod- 
ucts to be accumulated in the designated counters, the ma- 
trix can be divided into six (or 12) groups, with two (or 
one) elements of the B matrix in a card. On the other hand, 
if the elements are two-digit numbers, four elements of the 
final matrix can easily be obtained simultaneously. 



DISCUSSION 

Dr. Sherman: Could the described method of multiplying 
a vector by a matrix be extended to include multiplication 
of «th order matrices ? 

Miss Krawits: It could be extended very easily. It is only 
necessary to duplicate the number of cards in the B matrix 
as many times as there are rows in the A matrix. 

Dr. Sherman: Is there any simpler method that anyone 
has developed ? 

Miss Krawits: It is very difficult to answer your question 
directly. What may be the simplest method for one series 
of problems may be impossible for another. The simplest 
method, as I have stated before, depends upon the order of 
the matrix, the size of the elements, the machines available, 
and the number of multiplications in the problem. The 
method I have described is the most general, but is not the 
most efficient in all cases. For example, Dr. Petrie has used 
an entirely different approach. His problem was to multiply 
a 17th-order vector matrix by a 17th-order square matrix 
on the 602- A. All of the 17 elements of the vector are read 
into the machine on the first card cycle, and stored for the 
remainder of the calculation. Each element of the square 
matrix is put on a card (289 cards) and arranged in order 
of columns. After each element of a column of the square 



matrix is read by the machine it is multiplied by the corre- 
sponding element of the vector. The sum of the products of 
the vector by each column is accumulated and punched on 
a trailer card. This method is particularly efficient if a given 
matrix is to be used many times with different vector mat- 
rices. The fully equipped 602- A could handle a 24th-order 
matrix and vector, provided that the elements were four- 
digit numbers. 

Dr. Sherman: There is one slight thing that can be done 
to speed this process. It is not necessary, of course, to punch 
the a's into the card that contains the &'s. You can use one 
card as the master and the other as a detail card and, in 
general, eliminate the necessity for reproducing. 

Miss Krawits: It is true that the method you refer to 
eliminates the reproducing operation; it also keeps the origi- 
nal matrices intact, allowing them to be used over and over 
again. However, the machine time on the 602-A will be 
increased. The benefits of the two methods should be con- 
sidered with each problem. 

Dr. King: Have you ever modified this method to use 
complex numbers ? The principle seems to be applicable. 

Miss Krawits: Complex numbers can be handled just as 
easily as real numbers. It is true that the number of multi- 
plications and additions will be increased, but all of the 
methods described can be easily extended to include com- 
plex numbers. 

Dr. Hurd: One application of matrix by vector multipli- 
cation known to many of you is the application in spec- 
trometry. For a given instrument, you have a calibration 
which determines the coefficients in the left-hand side of an 
equation. With repeated routine analysis, the elements 
which change are the elements which are on the right-hand 
side of the equation. If the calibration is sufficient for three 
months and you have a good many routine analyses, it is 
profitable to invert the matrix and each time perform the 
matrix by vector multiplication. Similarly, many of you are 
interested in vibration analysis of mathematical models, 
which are in the form of linear ordinary differential equa- 
tions with constant coefficients. It is possible to write such 
an equation in matrix form. In order to obtain a set of 
fundamental solutions of the original equation, one way to 
approach it is to start with a vector and find out successive 
powers of the matrix itself times the vector, which again 
brings in the matrix by vector multiplication. In the air- 
craft industry where one problem is flutter analysis, this 
iterative method is used so much that some of our calcu- 
lators are busy 24 hours a day on this particular problem. 

Mr. Walker: In the problem where there are a great 
number of zeros in a matrix, if it were possible to arrange 
these zeros so that it approaches a triangular matrix, would 
that not greatly simplify the problem ? You could treat the 
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original triangular matrix in a method of elimination to 
take care of the few coefficients that do occur below the 
diagonal term and then invert the top triangle of the matrix. 

Dr. Grosch: Indeed, that is a fruitful idea. It is not al- 
ways easy to arrange the matrix in a purely triangular 
form, but it is very easy to use a matrix with a large number 
of zeros with standard equipment if you make a very careful 
sorting code. 

Dr. Sherman: In recalibrating from time to time there 
may be one critical calibration that has to be done often. 
You want to calculate the new inverse that results from 
just changing one column of the old one. I have worked 
with methods, and there are rather simple methods of cal- 



culating a new inverse from an old inverse as a result of 
making various changes, such as replacing one column by a 
new column or one row and column or one element, or 
deleting a row and column, or augmenting a row and col- 
umn. The practical situation is much more encouraging be- 
cause instead of having to go to the trouble of inverting the 
matrix, you can save a tremendous amount of work through 
a whole series of associated problems. We actually do that. 
In practice, you have a rather interesting situation of being 
able to take a large order inverse and make a rather large 
number of minor changes of rows and columns, or elements 
in rows, and maintain the original work of inverting that 
large matrix. 



Numerical Solution of Two Simultaneous 
Second-Order Differential Equations 

WALTER H. JOHNSON 
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THIS PAPER describes the numerical solution, using 
the IBM Type 604 Calculating Punch, of the following set 
of two simultaneous second-order linear differential equa- 
tions: 

e + J'ge-K^^ = 0, (1) 

i + <4<\>- Ked = Fit) . 



High 

Pressure 

Turbine 



Low 

Pressure 

Turbine 



AC Generator 







/0 






h 


ly 






e,ke 


<l>,k^ 







Figure 1 

These equations express the relations that exist between 
the angular displacements, acceleration of angular displace- 
ments, reaction torque, moments of inertia, and tortional 
stiffness of the shafts for a turbine-generator set, where: 

6 = angular displacement of the shaft between the 

high pressure and low pressure turbines. 
<f> = angular displacement of the shaft between the 
low pressure turbine and the AC generator. 

oyg — kff I — -}-y- ] = natural frequency of shaft 6. 
0)0 = ^0 ( 7- +7- ) = natural frequency of shaft <^. 

\i0 iy/ 



ke = torsional stiffness of shaft 6. 
k^ = torsional stiffness of shaft <f>. 

h 

h 

h 

reaction torq ue 



K, = ^ 



A-* = -r- 



F(0 = 



(Given for each 20°) 



Numerical Data 

50,000 kw AC generator set, tandem compound, double 
flow, 5-bearing. 
3,600 rpm 

h = 2,930 lb. in. sec.^ 
I^ = 9,5001b.in. sec.2 

/^= 15,650 lb. in. sec.2 (2) 

kg = 515 lb. in. 
k^ = 595 lb. in. 
f(t) is tabulated for each 20° of rotation from t — 0. 
20° 1 



^ ~ 360° 3600 
h = 20 degrees = 



60 = 0.000926 seconds 
20° 



57.2958°/radian 



(3) 



I = 2.86479/radian 
^ = 8.20702 



Approximation of Finite Differences 

F(t) is a known function of t over the range of integra- 
tion. Some constant increment in time, h, is decided upon 
so that 

tn + l = tn-\-h. (4) 

In this case h is taken to be the time required for the gen- 
erator to rotate 20°, since F{t) was given for this time 
increment, 

h = .000926 seconds. 
Suppose that in the region of an arbitrary point, f„, within 
the interval of integration, that the solutions 6 and ^ and 
their second derivatives can be developed into a Taylor's 
series. Then, 



1 + 1= Gn -\- hOn + :^ ^n 






+ 



4! 



+ 



and 



_ J h^" h^iJi h'^iv 



2r" ' 3!^" ' 4! 

II IT III }fi IV /j3 V /j4 VI 

also On + l = 6n -h Mn -{- Yi_ On + Y\ ^n -\- ^On -\- 

II II III /}2 IV /j3 V /i* VI 

and </>ra + l — <t>n -{- h(f>n -\- ^,4>n + 77</>n + TT </>« + 



(5) 



2! 



3! 



4! 
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A modified function can be defined which will eliminate 
all odd derivatives and the fourth derivative. This modified 
function is 

^ = ^-12^ (6) 

u 
Insert the expansions of ^n+i and Bn-^\ from equations (5) 

to obtain x^+i. Thus, 



^^i _ :^fl^ _ ^* fl^ _ :^fl _ 
12 " 12 " 24" 72" 



(7) 



or Xn+l = 

f'n + /^^„ + ^ e« + j2 ^« - 180 ^" ~ 480 ^" ~ 

By changing the sign of h, you can now write the expansion 
of the modified function at tn-i. That is; 



Xn—i = On ■— ndn + ^ 


U /j3 /// 


■^ 4! 


^n + ^ ^« 




// ;j3 /// 
^n + T2 ^n 


^4 
24 


JF /,5F 


or 

/ 5 


// 1 /// 


^180 


(8) 

A^V - 



A^',r„ =^n+i - 2;irn + Xn-l 



Similar series pertain for the modified function 3;. 

The central difference of the modified function x about 
the point tn is 

%_2 Xfi—2 
%— 1 Xn—\ 

A Xn — \/2 
tfi X<n, A Xfi 

^ •^» + l/2 

^n + 2 -*'ra + 2 
A^^^„ = A^^„ + i/2 - A^;r„_i/2 

^^ -^M + l Xn \Xn Xn,—\) 

(9) 



Substituting from (6), (7) and (8) into (9) gives 

1 



/ tr // 1 /// 



12 



12 



180 



h^On- ... 



n 



- 20^ + h^ On 



(10) 



/ C // 1 /// 1 F 

+ On- hOn +~h^On -y^k^ On +i^^' ^n 



or A^';ir„ = h^ 



1 ^^ 
240*°""- 



By neglecting the terms of the 6th order and higher, and 
thus forming a truncation error, we have 



likewise 



L'^Xn = h^ X 


-» 




^"yn = h^ 4>n 


_^ 


f A"y„ 



(11) 



By substitution into equations (6) from (11) we can state 
the values of the functions and </> in terms of the modified 
functions x and -y and their second differences: 



" "" " 12 V h^ )' 



and </)„ = y^ + 



Xn + A^^:r„ _ :r»-n + lO-Tn + ^n-1 

12 ~ 12 

12 12 



(12) 



Substituting in the original differential equations the ex- 
pressions derived for the functions and </> and their second 
derivatives in terms of the modified functions x and y, we 
obtain 

^^-l^yn + ^) - Ke (^Xn + ^") = j{tn) , 

(13) 



or 



^"Xn 



'"y\h+% 



Ka 



12 



= KeXn - iO%yn-\-f{tn) 



Solve for A^^;r„ and A^^3;„ to obtain the equations: 
A^^jf« = ^^„ + Byn+ Cj{tn) , 
A^% = £;r„ + Fyn + Gf(tn) , 



KffK^ — 12(0,9 



where A = 



V/t^^l2/ 



/ 



(14) 
B = hV 









F = 



G = 



KgK^ — 1^*^1(^+12) 
7 



J 



■' ~ ^ V^' 12/U' 12/ 12 ■ 
The coefficients A, B, C, B, F, and G are properties of 
the physical characteristics of the system only and remain 
constant for the entire integration — except for coefficients 
C and G, which contain the forcing function and must be 
computed for each integration step. 

Consider again the central difference formula: 

A"jir„ = Xn^l — 2Xn + Xn-l 

and A'^y„ = yn+i — 2y„ + y„_i , 

or Xn+l = A^^Xn + 2Xn — Xn-l (15) 

and yn+i = A^^y^ + 2y„ - y„_i . 
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Substituting into (15) from (14) for the second differences 
we have; 

^n + l = ^A'„ + Bjn + Cj{tn) + 2Xn - Xn-x , 
or Xn + l = (2-\-A) Xn + Byn + Cj{tn) - Xn-l 

and yn+i = BXn + fjn + Gf{tn) + 2yn - yn-1 , (16) 
or y„+i = (2+F) );„ + £.r„ + Gj{tn) - yn-i • 

Starting with initial values for Xq, x, yo, and y, the inte- 
gration can be carried forward to obtain the modified func- 
tions, x{t) and y{t). By employing equations (12), these 
solutions for the modified function are immediately trans- 
formed to the functions B{t) and ^{t)\ 



i.e., Bn = 



12 



(12) 



The constants for the problem are 

{2 + A) = -i-0.59410, 
B = -f 0.36467 , 
(2-f F) = H- 1.34706, 
B = -f 0.31564 . 
The coefficients Cf{tn) and Gf{tn) are tabulated for inte- 
gration points tn. 

Solution on the IBM Type 604 Calculating Punch 

Cards are prepared in pairs; i.e., one card each for the 
$ and ^ shaft for each point tn (Figure 2). 

Card Form 



col. 1-2 



Cardl 
Card 2 



col. 3-7 
BorE 



.36467 
.31564 



col. 8-13 

(2-f J) 

or (2-hF) 



.59410 
1.34706 
etc. 



col. 14-18 

cm 

or Gf{t) 



.00015 
.00481 



col. 19-23 



col. 24-28 



Figure 2 

Assignment of Data in 604 Calculator 

col. 3-7 = B or B assigned to MQ (Multiplier Quotient 

register) 
col. 8-13 = (2-]- A) or {2-j-P) is assigned to FS 1 and 3 

(Factor Storage units) 
col. 14-18 = Cf{t) or Gf{t) is assigned to GS 1 and 3 
(General Storage units) 



Sequence card ^„: 

1. RO GS 4, multiply -I- 

2. ROGS2, RIMQ 

3. RO FS 1 and 3, multiply -f 

4. RO GS 1 and 3, RI EC + 6th 

(2 + A)Xn-{-Byn+Cf{tn) 

5. RO FS 2, RI EC— 6th {2 -\- A)xn + By^, 

+ Cf{tn) — Xn-l ~ ^n + 1 



Byn 

Xfi 

{2 + A)xn-\-Byn 
Cfitn); 



6. ^ adjust 5th 

7. RO and RS EC 6th, RIMQ 

8. Emit 1 + 2nd 

9. Emit 2 + 1st 

10. RO and RS EC, RI GS 1 and 3 

11. ROMQ,RlEC+ 1st 

12. RO GS 2, RI EC + 2nd 

13. RO FS 2, RI EC -f 1st 

14. ROGS2, RIFS2 

15. ROMQ, RIGS2 

16. RO GS 1 and 3, divide 6n = 



■*■» + ! 

12 
12 

■^n+l 

•^n + 1 + 10x„ 

Xn+i + 10;i:n + Xn-i 

shift Xn to ;r„_i 
shift -rw+i to j*r„ 

■r„+l — 10.yn — Xn-l 

12 



DISCUSSION 

Mr. Brown: I would like to ask about the possibility of 
getting these difference equations in a more direct manner. 
What is done here is to start with a differential equation and 
then go through an elaborate process of making approxi- 
mations in terms of differences. Then, finally, one finishes 
with a set of difference equations. Now, in this case, the 
original differential equations themselves can be obtained 
by varying an integral. Instead of going through all the 
manipulation with the differential equation, one could con- 
ceivably start with an approximation in terms of differ- 
ences for the kinetic and potential energy terms in Hamil- 
ton's integral. Then, carry out a variational process and 
obtain difference equations without ever handling the dif- 
ferential equations and without having to put all these ap- 
proximations into the differential equations themselves. In 
many physical problems which can be expressed as the 
problem of varying an integral, the usual procedure is to do 
something like that; for example, to introduce a polynomial 
or other approximation with undetermined parameters in 
the integral and then to vary the parameters to minimize 
the integral. That would be the logical thing to do here. Has 
that been looked into in this type of problem, and if so what 
are the relative merits of each method ? 

Mr. Johnson: We did not look into that. I presume you 
could do it as you have outlined. 

Mr. Collins: It is unnecessary to solve this problem by 
approximate methods as far as I can tell. If you have con- 
stant coefficients, it looks as though you could obtain two 
quadratic algebraic equations and reduce them to a fourth 
order polynomial and calculate the necessary roots. You 
would have to write out an explicit formula for the answer, 
and then evaluate all the points you wanted to from trig- 
onometric exponential functions. I wondered why it is de- 
sirable to do the problem by this approximate method. 
There probably is a good reason. 

Mr. Sheldon: I think it is really easier to solve this prob- 
lem by the integration of the differential equation instead of 
evaluating the integral you would obtain if you solved it 
explicitly. 



Numerical Evaluation of Integrals of the Form j f{x)g{x)dx 
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SUPPOSE we wish to evaluate numerically 

j{x)g{x)dx , 

and suppose g{x) varies much more rapidly than j{x) so 
that much finer intervals are required to represent g{x) 
by an interpolating polynomial than would be required for 
j{x). It is then worth while to consider calculating special 
Lagrangian integration coefficients, di, which take account 
of the variation oi g{x) so that 

j(x)g(x)dx = y dijt 



£ 



with the ordinates fi chosen at intervals in x which need 
only be fine enough to represent f(x) by an interpolating 
polynomial. Cases where this approach may be especially 
useful are: 

1. f(x) depends on one or more parameters, and inte- 
grations are desired over ranges of the parameters. Then 
the number of operations required to solve the problem may 
be very substantially reduced by this technique. 

2. An integrand, P{x), has a singularity g{x) where 
g{x) can be integrated analytically, and jg(x)dx is finite. 

Then Fix) =^9i^) = fi^)9i^) • 

3. g(x) is periodic and varies so rapidly that we can 
represent P{x) by an interpolating polynomial with inter- 
vals in X of one period of g(x). In this case we obtain 
especially simple integration formulae. 

Derivation of Formulae 

for the Special Integration Coefficients 

For simplicity, let the range of integration be to a. Let 
(0,a) be divided into n equal intervals (0,h), (h,2h), . . . , 
[(n-l)h,nh].Then 



/f{x)g(x)dx = I 



h/2 



f{x)g(x)dx-\- 



n—l /^U+l/2)h /^nh 

^ / f{x)g{x)dx\ / f{x)g{x)dx 

}=1 •/(i-l/2)A t/{n-l/2)h 



(1) 



Let X = h{j-\-$) in the /th interval. 



Then 



/o / X»l/2 

f{x)g{x)dx = h\ j fihi)gihi)di 



n—X /^l/2 
j=lt/-l/2 

+lf[h(n-\-$)]g[h{n+i)]dA 
J-l/2 )(2) 

TO 

Let* f[h{j-j-i)] = \ bkji^' where w is a positive even 
number, A (3) 

TO/2 

and where &&,- = > aikfi+i. (4) 



?==— TO/2 



Then,/ f(x)g(x)dx 




'It''' f 

\k=0 »/o 

m /^ 

fc==0 •/—I, 



1/2 

$''g[h$]di (5) 



n — 1 TO /»l/2 TO /^O 

+ > > leg[h(j+i)]di+}bjcnieg[h(n+i)]di}. 



Let Cofc= / e9{hi)di 



I 



1/2 



aHere is a representation in Lagrangian form of an interpolating 
polynomial. I have found that in most cases the Lagrangian form is 
most convenient. An exception is case 3 above, where it is best to 
introduce interpolating polynomials with coefficients given in terms 
of differences. The particular interpolating polynomial introduced 
here involves ordinates outside the range of integration. Of course, 
if f(x) does not exist outside the range of integration, different 
interpolating polynomials must be used near the ends of the range 
of integration. 
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Cjk 



= leg 

J-X/2 
•/-1/2 



9[Ki-\-md^ 1= 1,2, 



n— 1 
(6) 



1/2 

n w m/2 



/a n m m/2 

f(x)g(x)dx = h\ \ y aikCj,cfj+i 
}=0 k=0 l= — m/2 



(7) 



/ 



Let / = />-/. Then, 

m/2 p+m/2^n m 

j{x)g{x)dx = ^ 2j X X ^^-^•^^^^^'' 



p=— TO/2 j=p-m/2^0 k^\p—j\ (8) 

p+m/2^n m 



SO that 
and 



dp = h y y ap—j^kCjjc , 



/=sp-TO/2^0 fe^lp— il 
p=m+m/2 



/o p=Jt-fTO/2 

f(x)g{x)dx = y rfj^ 
p=— OT/2 



(9) 



(10) 



After the quantities dp have been evaluated, the integral 
can be evaluated efficiently by machine for any choice of the 
function f(x), provided only that f(x) can be represented 
by interpolating polynomials in the interval chosen. The 
evaluation of the coefficients dp is often the central problem. 
First of all one must obtain the coefficients Cjk. Sometimes 
the Cjjc can be obtained by analytical integration, in which 
case they are represented by formulae which may be evalu- 
ated by hand or by machine. Otherwise, the Cjk must be 
obtained by numerical integration. Once the Cjk are found, 
the dp's may be evaluated from equation (9). This summa- 
tion is somewhat tedious to carry out if the card volume is 
small. Equation (9) is a formula for the dp's which gives 
results with accuracy analagous to Bessel's integration 
formula in the ordinary theory of numerical integration, 
and it is the most accurate type of formula we can obtain 
using equal intervals and interpolating polynomials. A 
simpler procedure which we have used, but which is not so 
accurate, involves strip formulae and Vandermonde mat- 
rices as follows: 

Divide the range of integration into subintervals (0,h), 
(h,2h), . . . , [(n—l)h,nh] as before, but now let 



/a n /m— 1 /^m(p + l)h 

f(x)g(x)dx =2.1 f{x)g{x)dx 
p=0 %/mph 



(11) 



Vntx = h{m{p-\-\/2) -|- |] so that 



/a n/m—1 /*i 

f{x)g{x)dx = h\ I j 

/ 



n /m— 1 •» m/2 

f{x)gix)dx = hY I f{h[m(p-]-l/2) + i] } 

-m/2 

' g{h[m(p+l/2) + ^]}d^ , (12) 

• m/2 

and let I f{h[m(p+l/2) + $]} 

'—m/2 

j=m/2 

' g{h[m(p-\-l/2) + $]}d^ = V dmp+jfmp+j. (13) 

^=— TO/2 

In order to evaluate the coefficients dmp+j, let 
f{h[m{p -{-1/2) + ^]} equal successively 1, ^, 1^ . . . , |^ 

••to/2 

and let / eg{h[m{p+l/2) + i]}d$ = dp. (14) 

»/— m/2 

Then we obtain the system of (w+1) simultaneous equa- 
tions 

y=TO/2 

c^ = > ^jdmp+j w = 0,1,2, . . . , w (15) 
i==— TO/2 
which can be solved for the dmp+j, obtaining 

ra=OT 
dmp+J= V(I7)-^Ctop. (16) 

M=0 

After we have calculated d^p+j from (16) we must add 
together rfTOp+ TO/2 and dm(p+i)-m/2 as they are both coeffi- 
cients of the same ordinate. Here we have developed "strip" 
and "repeated strip" formulae which are analogous in ac- 
curacy to such formulae as Simpson's rule and Weddle's 
rule in the ordinary methods of numerical integration. 

Applications 

1. Evaluate / f(x)cos kx dx when k is so large that 



/ 



we may take steps of one period in the cosine and still 
represent f{x) accurately by an interpolating polynomial. 
In this case, we use the more accurate formula (9) for the 
di's and Stirling's interpolation formula in terms of differ- 
ences in place of equation (3). 

Because the average value of the cosine is zero, all the 
"zero-order" terms in the interpolating polynomial con- 
tribute nothing to the integral. Because the cosine is an 
even function, all the odd order terms cancel. We are left 
with a sum over the even order terms which are even order 
differences of f(x) times constant coefficients. These are 
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integrated to odd-order differences at the beginning and at 
the end of the range of integration. Taking account of all 
differences up to and including the 6th, we obtain 



f* CO 

j{x) COS kx dx = U 0.0016007662 ["/ (=^^ - /(y)1 



(17) 




-0.0186667286 f/ (-^\ - /(y)1 

+0.2071512979[/(:^)-/Q]}. 
For example, evaluate 



/ 



^~^cos SOttx dx . 



We obtain from (17) 

g-^cos SOTT.r dx = 0.0000405268311. 
By analytical evaluation we obtain 

/oo 
^-^cos 50^;ir dx = 0.0000405268310. 

To obtain this result using, for instance, Simpson's rule 
would require us to add many thousand values of the 
integrand. 



2. Evaluate 



/ 



g-a)/2 Ko(:r/2)cos kx dx where Ko{x) 



is a Bessel function of pure imaginary argument, defined in 
Watson's Bessel Functions} Near the origin Ko{x/2) ;^ 
— y — loge(.s;/4) where y is the Euler-Maclaurin constant. 
Dividing by logglO, we have 
Ko{x/2) y 



Let f{x) 



logelO logelO 

e~'^''^Ko(x/2)cos kx 



- logio(^/4) 



0.48. 



-y/logelO - logio(^/4) > for ^ _ 
g{x) = -y/logelO - logio(^/4) ) 

Then f(x) is "slowly-varying" as compared to p(x). To 
evaluate the integral, take steps of .04 in x and use m = 4 
(five-ordinate strips) in equations (14), (15) and (16). 
Then, 



C4p 



=/ 



^{-y/logelO - logio[.16(/'+l/2) + .04i]]d$ 



For p — we have 

cO = .2530485 
cl = .0694871- 
cl ^ .3682811 
cl = .0926495- 



The matrix (0) is 










/-I +] 

-2 -] 


L 1 ] 
L ] 


L 1\ 
L 2 






+4 ] 
-8 -] 


L ] 
[ ] 


L 4 
[ 8 




\ 16 ] 


L ] 


L 16^ 


> 


and we find that (0^ 


-1 is 








/o 



1 



\o 


1/12 

-2/3 

2/3 
-1/12 


-1/24 
2/3 

-5/4 
2/3 

-1/24 


-1/1^ 
1/6 

-1/6 


I l/24\ 
-1/6 
1/4 
-1/6 

I 1/24/, 



so that from formula (16) we obtain 
rf_2 = .0252662 
(/_i= .1216794 
rfo = .0247837 
dr = .0599131 
d2 = .0214058 . 

Repeating this process for p — 1,2,3 we obtain 

rf_2 = .0252662 

rf_i = .1216794 

do = .0247837 

rfi = .0599131 
-^2 = .0355995 

dz = .0600419 

^4 = .0203343 

^5 = .0516705 
-rfe = .0209827 

di = .0453289 

ds = .0158488 

^9 = .0403588 
-dio = .0083192 . 
idmp+m/2 and dm{p+i)-m/2 have been added together.) 

For the range .48 to oo, which we take to be .48 to 10.0, 
we use Gregory's interpolation formula in Lagrangian form, 
accurate to 6th differences, and continue integrating in 
steps in x of .04. We use the integrand e~''^^Ko{x/2)cos kx 
without splitting it into product functions. 

The integral can be evaluated analytically.^ 



/ 



"/^i^o i^/2) cos kx dx— 2 • real part 



) cos~^(l— 2i^) 
\\/l-{\-2izy 



We obtained the following results:^ 



d = .9283464 



bThe evaluation of this integral was undertaken in connection with 
a problem at the IBM Technical Computing Bureau in New York, 
where we were really interested in the values of a somewhat simi- 
lar integral which could not be evaluated analytically. The above 
integral was evaluated as a check on the theory and numerical work 
in obtaining the d's, which were then used for the actual integral 
which interested us. 
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k 


By Numerical Integration 


Correct Value 


0.0 


\ 99997 S 


2.00000 


0.1 


1.989402 


1.98963 


0.2 


1.958604 


1.95858 


0.4 


1.847696 


1.84764 


0.7 


1.616215 


1.61902 


1.0 


1.380436 


1.38051 


2.0 


0.843857 


0.84383 


4.0 


0.437586 


0.43749 


7.0 


0.246750 


0.24624 



For larger values of k, we have developed special inte- 
gration coefficients for the function 

g{x) = [y/loge 10 - logio (^/4)] cos kx . 
For this choice oi g{x) the c's themselves had to be evalu- 
ated by numerical integration. 

In the application of this method to various problems, I 
acknowledge with pleasure the assistance of Mr. Liston 
Tatem and Mr. Daniel Ladd of the Applied Science Depart- 
ment, International Business Machines Corporation. 

rijfrrence; 

1. G. N. Watson, Theory of Bessel Functions, Second Edition 
(MacMillan Company, 1948), pp. 388. 

DISCUSSION 

Mr. Bejarano: Can this process be used for triple inte- 
grals ? 



Mr. Sheldon: Yes. 

Dr. Buchhoh: Is it still profitable to use this method if 
you must resort to numerical integration to obtain the c 
coefficients, instead of using analytical methods ? 

Mr. Sheldon: Yes, that is what we are doing in the case 
where g{x) = cos kx log x dx. The function can be inte- 
grated analytically, but it involves the cosine integral which 
is a tabulated function. Rather than evaluate the c's analyti- 
cally, in this case it is simplier to take very fine intervals 
and evaluate the CjkS numerically. 

Mr. Clamons: I have just one comment to make on this 
subject. Many times you run into a problem, in which there 
is a product function like this, in connection with the ex- 
perimental data where one function is experimental and the 
other is analytical. It pays to establish the analytical func- 
tion by means of that integral keeping the g function on the 
inside and then establishing a matrix deck. 

Dr. Hurd: I think this whole discussion illustrates again 
the fact which we know very well, namely, that computing 
machines do not replace mathematicians. Here is an in- 
stance in which for a given integral, such as Mr. Sheldon 
has discussed, there are obvious methods of treating them 
directly by the definition of the integral. If you do this, you 
have a very slow computing process and the computing 
machine turns for days. If you apply some mathematical 
knowledge in advance, you cut the amount of computation 
time down. I am always glad to think that computing ma- 
chines make the work of mathematicians more valuable. 



The Use of Orthogonal Polynomials in Curve 
Fitting and degression Analysis'^ 
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IN THE PROBLEM of two-variable curve fitting, a 
simple polynomial in powers of the independent variable is 
undoubtedly the most frequently used function in cases 
where the form is unknown from empirical or theoretical 
knowledge of the problem. Denoting the independent vari- 
able by X and the dependent variable by Y , the polynomial 
may be written as 

Y = ao + aiX -H 0-2X2 + .... (1) 

The parameters a^, ax, a^, .... are to be evaluated by the 
method of least squares; i.e., so that the sum of the squares 
of deviations between the calculated and observed values 
of F is a minimum. 

For the case that X values are without error and the Y 
values have equal weightings (precision), the so-called 
normal equations for the evaluation of the parameters are 

TVao + (2X )ax -f {■%X^)a^ +.... = 2F (2) 

{%X )flo -f- (2X2)ai + {■%X^)a2 -f- . . . . = ^XY (3) 

(5X2)ao -h {%X^)ax + (5X4)a2 -f- . . . . = 5X^7 (4) 



In the above equations, the index and limits of summa- 
tion are omitted, for in each case the summation is carried 
out over all the values of X and Y — i.e., from 1 to N. 
A more general treatment of the least squares problem, in 
which both the dependent and the independent variables 
may be subject to error and in which the parameters may 
enter the equations nonlinearly, has been published.-^ The 
matrix of the coefficients of the parameters in a set of nor- 
mal equations is always symmetrical about the main di- 
agonal. These equations may be solved by any of the stand- 
ard methods of solving linear simultaneous, algebraic 
equations. 

The well-known method of curve fitting outlined in the 
foregoing paragraphs has two undesirable characteristics: 

*This paper was presented informally. 



1. The computational labor of solving the normal equa- 
tions becomes considerable when the number of 
parameters is large (e.g., greater than 5). 

2. If the polynomial of the ^th degree is fitted by the 
method of least squares and it is then decided to add 
the additional term aitj^iX^^'^, not only must a^+i be 
evaluated but all the parameters as well, because their 
values will change by the inclusion of this additional 
term. Thus, in the usual method of curve fitting, the 
degree of the polynomial must be decided at the 
outset. 



Orthogonal Polynomial Method^-^ 

Instead of expressing F as a polynomial directly in 
powers of X, it may be written more generally as 

F = ^0 + Axti + ^2^2 + ...., (5) 

in which ^1 denotes a linear function of X, ^2 a quadratic 
function of X, etc. 

If the I's are so chosen that 



I 



UXic) li(x,) = 



(6) 



for all values of i and /, except when i = /, the polynomials 
are termed orthogonal. The usefulness of this property of 
orthogonality in curve fitting lies in the fact that all the 
coefficients of the parameters in the normal equations, ex- 
cept those on the main diagonal, are zero; thus, a complete 
separation of the parameters is achieved. 

The normal equations then become: 

NAq = 2F (7) 

(S|f)^x = 2^iF (8) 

{Ml)A2 = ^^-2Y (9) 
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For the case that the values of X are equally spaced, the 
orthogonal polynomials can be easily derived.* It is this case 
which will be discussed below. 

In practice, it is more convenient to deal with the equa- 
tion in the form 

F = ^(, + ^ili + ^^|^ + .... , (10) 

in which the ^' polynomials are related to the | polynomials 
by means of equation 

^i^Ui. (11) 

The \i values are so chosen that the |' values are integers 
reduced to the lowest terms. The relationships between |' 
and X up to the fifth degree are: 



X 



Y 



li 



|2 



ri = At {X-X) , 



I2 = A2 Y 



(x-xy 



N^- 



12 



^1- 



(12) 
(13) 



iz = Xs[(X-Xr - (X-Z)(^^^)],(14) 






. 3(JV^-l)(Ar^-9) -| ^j3^ 



560 



^5 = A5[(X-Z)« - (Z-Z)3^^^^^^-^^ 

}< 



+ (X-X)-^^^ 



18 
230N^ + 407 



These equations enable one to transform the equation in 
the i's to an equation directly in terms of X, if so desired. 
However, the sum of the squares of the deviations from the 
regression equation may be obtained directly in terms of ^'s. 

The equation is 

2A2F = 2F2 - A'o^(Y) 

-A'a{Yi[) -Aa(Y$'2) -... . 

The first term, SF^, is the sum of the squares of the F 
values about zero. The second term, ^o5( F), is the amount 
by which the sum of squares about zero is decreased when 
taken about the mean. The third term, ^iS( F|i), gives the 
amount by which the sum of the squares is further decreased 
when the deviations are taken about the best (least squares) 
straight line, etc. 

Illustrative Numerical Example 

As a simple numerical example illustrating the applica- 
tion of orthogonal polynomials to the determination of a 
regression equation, consider the X and F values given in 
the first two columns of the table which follows : 



4 


1.82 


-3 


5 


5 


6.13 


-2 





6 


12.09 


-1 


-3 


7 


19.47 





-4 


8 


29.80 


1 


-3 


9 


42.12 


2 





10 


55.91 


3 


5 



28 



84 



5^? = 
If the A' and F values are fitted to a polynomial 

Y = A + BX+CX^, (18) 

by the standard method of least squares the normal equa- 
tions are: 

7A+ 495 -f 27C = 167.34 (19) 

49^ + 371B + 2989C = 1423.34 (20) 

371^ 4- 29895 + 25235C = 12481.56 . (21) 

The solution of these equations yields the polynomial 

F = 6.50494 - 5.18474Z + 1.01309^^ . (22) 

To obtain this same result by the application of or- 
thogonal polynomials, the following quantities have been 
calculated: 



% = 23.90571 

N 



(23) 

2(^{F) = 251.96 (24) 

Mi2Y)= 85.1 (25) 

Dividing the latter two steps by the number standing at the 
bottom of the corresponding $' columns, the following 
parameters are obtained: 



^; = 2^=251^ = 8.99857 



A^ = 



M^V) 



85.1 
84 



= 1.01309 



(26) 



(27) 



The corresponding orthogonal polynomial expression for 
Fis 

Y =^ 7 + A[ii + A'^i^ (28) 

= 23.90571 + 8.899857 H + 1.01309 ^'2 . 

To transform the above equation in $' into an equation in X, 
the table of relationships previously given is used. (In this 
example Ai and A2 are equal to 1, so that |i = |i and 
^2 = $2- Also, 7 is the mean value of X.) 

Y = 23.90571 + 8.99857(X-7) 

+ 1.01309 [(Z-7)2-4] (29) 
= 19.85333 -f 8.99857 (X- 7) 

+ 1.01309 (Z-7)2. 
The above equation in powers of X — 7 may be converted to 
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an equation in powers of A^ by a method employing syn- 
thetic division. The result is 

Y = 1.01309X2 - 5.18476X + 6.50499 . (30) 

A comparison of equations (22) and (30) shows that 
they are identical to within a few units in the last figure. 

In order to obtain the calculated values of Y or to carry 
out a regression or correlation analysis, it is not necessary 
to convert the equation in |' to an equation in X. The details 
need not be given here.^ 

Orthogonal polynomials may be used very efficiently to 
obtain an wth degree polynomial corresponding to a tabu- 
lation oi N—1 equally-spaced values of F. 
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DISCUSSION 

Mr. Brown: As Dr. Sherman remarked, these tabular 
polynomials exist in convenient form only when the points 
are equally spaced. Experimentally, and on a small scale, I 
have tried calculating them and using them in a few cases 
when the points were not equally spaced. In that case, of 
course, the usual procedure is to form the normal equations 



of B'-^ and go ahead and solve them. I have tried this alter- 
native method which consists in calculating, for the particu- 
lar set of given points, as many of the orthogonal polyno- 
mials as are needed for the particular problem. I am not 
sure about the relative amount of work involved in the 
conventional procedure, and in this procedure, but I am 
favorably impressed by what seemed at first glance a re- 
markably laborious task, namely, calculating a whole set of 
functions before you start your problem. I think the amount 
of work is about equal by the two methods. 

The method has a great advantage. If you want to fit 
more than one function to the same set of points, you may 
perhaps have eight pages of work for the first function, but 
seven of those pages consist in calculating orthogonal poly- 
nomials, and the eighth consists of evaluating the coefficients. 
If I now want to fit a second function, there is one more 
page, not eight more. If it were calculated by the conven- 
tional method, that would be a great deal of calculation, but, 
on the whole, I am impressed with that way of doing it. 
I think it merits further investigation. I don't have any final 
evaluation of the two methods. 

In some cases two variables may have almost the same 
effect on the variable, and they vary closely with each other, 
in which case there will be difficulty in solving because of 
. an almost vanishing determinant. The two variables are so 
closely correlated that there may be some indeterminateness 
in determining these coefficients. On the other hand, certain 
variables may have so little effect on the variable Y that 
they can probably be left out. I think it would be interest- 
ing to investigate the usefulness of orthogonalizing such 
functions. 



General Purpose Ten-Digit Arithmetic on the 
IBM Card-Programmed Electronic Calculator"^ 
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THE IBM Card-Programmed Electronic Calculator 
makes available to the computation laboratory of moderate 
size a computer of considerable flexibility. By means of ap- 
propriately designed control panels, it is possible to exploit 
the facilities for control provided by the accounting ma- 
chine, and thus to employ the equipment for the efficient 
solution of a wide variety of problems. 

We may conceive of two diametrically opposed methods 
for the utilization of the card-programmed electronic cal- 
culator. The control panels may be wired in such a manner 
as to permit the computer to perform the elementary arith- 
metic operations under control by suitable code punching in 
program cards which are read by the accounting machine. 
A given computational routine is then performed by reduc- 
ing the problem to a series of arithmetical operations to be 
performed in sequence, and the program for the problem 
consists of the sequence of coded instructions punched into 
cards, by means of which the computer is given step-by-step 
instructions for carrying out the desired routine. When em- 
ployed in this manner, the computer may be said to be a 
general-purpose computer. On the other hand, it is fre- 
quently possible to control simple sub-routines by means of 
control panel wiring, and it may be possible to generate 
some or all of the program within the machine instead of 
employing card control. In the extreme case, the computer 
may then be a single-purpose computer and may be said to 
be controlled by internal program. 

Substantial savings of operating time and greatly in- 
creased efficiency of operation frequently result from the 
specialization of the control circuits to a particular problem. 
In this way the substantial facilities for control and storage 
provided by the accounting machine and the subsidiary 
storage unit can be employed in conjunction with the elec- 
tronic calculator for the rapid solution of a particular prob- 
lem. However, a considerable amount of time is normally 
required for the design and installation of the required con- 
trol panels. This type of operation is, therefore, not usually 
justified unless the problem is of considerable magnitude. 

*This paper was presented informally. The supplementary notes 
were added subsequent to the presentation. 



Although the advantages to be gained by specialization 
of control panels should be examined for any particular 
problem, it will usually be the case that the average compu- 
tation of moderate duration can be most expediently carried 
out on a general-purpose computer. In this way, the design 
and installation of control panels for the problem at hand is 
entirely eliminated, and the planning for the problem con- 
sists only of the formulation of the card program routine. 

In view of the fact that the general-purpose computer is 
to be used in a wide variety of problems, it is desirable that 
it offer to the programmer the greatest possible flexibility 
consistent with the logical design of the equipment. It is 
desirable that it be possible to carry out the arithmetical 
operations with the minimum number of program cards. 
Furthermore, it is desired that facilities be provided for the 
automatic selection of alternative computational routines 
and for the selection of input data. 

The electronic storage capacity of the IBM Type 604 
Electronic Calculating Punch is assigned to provide four 
storage groups for ten-digit numbers. The units connected 
to channel A are designated electronic storage A. Similarly, 
electronic storage B and electronic storage C consist of the 
units connected to channel B and channel C, respectively. 
Auxiliary temporary storage is provided by additional units 
and is designated electronic storage D. 

Provision is made for the eight arithmetic operations 
summarized in Table I. These consist of the elementary 
operations, square root, and certain combinations. 

Table; I. Arithme;tic Orde;rs 





Order 


- 




Storage 


Result 




Code 


A 


B 


C 


D 


Decimal Point 





A = C 




C 





C 


D 


c =^ a 


1 


A + B = C 




c 


B 


c 


D 


c = a=:b 


2 


A - B = C 




c 


B 


c 


D 


c ^= a:= b 


3 


A X B = C 




c 


B 


c 


D 


c = a -}- b 


4 


A -^ B = C 




c 


B 


c 


A 


c = a — b-\-l 


5 


VA = C 




c 


C/10 


c 


A 


c = a/2 


6 


A + B + D = 


C 


c 


B 


c 


D 


dz=z c = a — b 


7 


A - B + D = 


C 


c 


B 


c 


D 


d = c — a= b 


8 


A X B + D = 


c 


c 


B 


c 


D 


d = c = a 4- & 
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81 
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The table gives the condition of each of the electronic stor- 
age groups at the conclusion of the operation. The factors 
are identified by the label of the electronic storage group in 
which they are stored or by the channel over which they are 
transmitted. The decimal point rule for each operation is 
also given, where a denotes the number of decades between 
the leftmost decade of the number A and the decimal point, 
and b, c, and d are similarly defined for the numbers B, C, 
and D, respectively. The order code zero results in a trans- 
fer of the factor A to electronic storage C. 

The result C of every operation is transferred to elec- 
tronic storage A, where it may be employed as a factor in 
the succeeding operation. The factor B is preserved by every 
operation and may be employed in succeeding operations 
unless the square root operation intervenes. The number in 
electronic storage D is unaffected by all operations except 
division and square root, in which cases it is replaced by 
the factor A. 

The provision for the retention through a sequence of 
operations of significant factors within the electronic stor- 
age of the computer contributes substantially to the effi- 
ciency and flexibility of the computer. It is frequently pos- 
sible to carry out routines of considerable length before 
sending results to external storage. These provisions effec- 
tively increase the over-all storage capacity of the computer, 
a consideration that may be significant in problems for 
which storage capacity is at a premium. 

The square root operation utilizes an iterative sub-routine, 
programmed on the 604 control panel. A special method of 
obtaining the first trial value of the root has been employed 
so that the complete operation requires a maximum time of 
two machine cycles. 

In addition to the transfer of the result C to electronic 
storage A for every operation, provision has been made for 
three optional transfers between the different electronic 
storage groups. These transfer orders follow the arithmetic 
operations. They control access to electronic storage D and 
supplement the provisions of the arithmetic operations for 
the destination of the operation result. The optional transfer 
orders are summarized in Table II. Except as noted, they 
can accompany any of the arithmetic orders. 

Tabi^e II. Trans:pe;r Orders 









Storage Result 






Code 


Order 


A 


B 


C 


D 


Notes 


6 

7 
8 


C->B 
C-^D 


C 

C 
C 


D 

C 
B 


C 

C 
C 


D 

D 
C 


not with arith. orders 

4, S, 6, 7, 8 
not with arith. order 5 
not with arith. order S 



The card-programmed calculator is regularly equipped 
with a shift unit that is associated with channel C. This unit 
makes it possible to shift a result C up to five decades to the 
left. The ease with which various problems can be pro- 
grammed is markedly increased if provisions are made for 
shifting on channel A of the factor A before it enters the 
computing unit. The field selector of the accounting machine 
can be employed to provide a shift to the left of up to four 
decades, or a shift to the right of up to five decades. The 
shift orders are summarized in Table III. 

Tabi^e; III. Shift Orders 



Code 
Channel A 


Result 


Code 
Channel C 


Result 





A 





C 


1 


\0A 


1 


IOC 


2 


WA 


2 


lO^C 


3 


WA 


3 


lO^C 


4 


K^A 


4 


10*C 


5 


10-M 


5 


lO^C 


6 


10-M 






7 


10-3^ 






8 


10-M 






9 


10-M 







The optional transfer orders contribute significantly to the 
flexibility and efficiency of the computer. 



The storage addresses of factors A and B and the address 
of the storage unit to which the result C is to be sent, to- 
gether with operation and transfer codes and shifting in- 
structions, constitute a ten-digit order code or word that 
normally is read by the accounting machine from columns 
4 to 13 of an order card. The location of the various parts 
of the complete order is shown in Figure 1. Input to the 
computer is achieved either by entry from the order card 
into channels A and B under control of a 00 address for the 
channel involved, or by an order termed spread load (SL) 
under X control, permitting direct entry from the order 
card into six of the seven accounting machine counter 
groups. The card fields associated with each counter group 
during a spread load operation are indicated in Figure 1. 
The punch unit is wired to punch from certain accounting 
machine counter and electronic storage groups as shown in 
Figure 1. In both input and output, negative numbers are 
punched as true figures with an X in the units position of 
the card field, and sign control within the computer is 
automatic. 

The provisions that have been described for the arith- 
metic, shift, and transfer operations are sufficient to permit 
the construction of efficient programs for computational 
routines involving purely algebraic operations. However, 
unless further implemented by provisions for program con- 
trol, it is necessary that the machine operator perform all of 
the discriminatory operations involved in the solution of the 
problem at hand. In many problems — as, for example, those 
of iterative nature — it is necessary to provide alternative 
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INSTRUCTIONS 



CHANNEL A 
ENTRY 



CHANNEL B 
ENTRY 



INSTRUCTIONS 



ORDER CODE and CHANNEL ENTRY 



COUNTER C 



COUNTER 2 



ooff 

1 2 3 
1 1 1 

222 
333 
444 
555 
666 
777 



999 

1 2 3 



0000000000 

14 15 16 1718 19 20 212223 
1111111111 

2222222222 
3333333333 
4444444444 
5555555555 
8666666666 
7777777777 



9999999999 

14 15 16 17 18 19 20 21 22 23 



0000000000 

2425 28 27 2829 30 313233 
1111111111 

2222222222 
3333333333 
4444444444 
5555555555 
6666 666666 
7777777777 



9999999999 

2425 26 27 282930 3132 33 



99 

34 35 



00 

37 38 
1 1 

22 
33 
44 
55 
66 
77 



99 

37 38 



00 

3940 
1 1 

22 
33 
44 
55 
66 
7 7 



99 

3940 



0000000000000000000000000000000000000000 

41 42 43 44 45 46 47 48 49 so 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 78 77 78 79 80 
1111111111111111111111111111111111111111 

222222222222222222222222222222222222222 2 
33333333333333 33 33333333333333333333333 3 
4444444444444444444444444444444444444444 
5555555555 5555 5555555555555555 555 5555 55 5 
66666666666666 6 6666 666666606666666666666 
77777777777777 7 7777 777 77777 777 7777777 77 7 
8888888886888888 8888888888888888888888 8 8 
9999999999999999999999999999999999999999 

41 42 43 44 45 46 47 48 49 50 51 S2 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 76 79 80 



Figure 1. Instruction Card 



computational routines, and at some point in the calculation 
to select the appropriate routine on the basis of criteria de- 
veloped by the calculation itself. 

A considerable amount of automatic program control can 
be achieved by the addition of two program control orders. 
The ten-digit instruction word is normally contained in col- 
umns 4 to 13 of the order card, and this field is designated 
the normal instruction field. An alternative field, called the 
transfer instruction field, is located ( Figure 1 ) in columns 
34 to 43 of the order card. If, at some point in a calculation, 
the number occupying electronic storage C is negative, and 
if "conditional transfer" (CT) is ordered at this point by 
an X punch in the order card, the computer will transfer 
from the normal instruction field to the transfer instruction 
field, and it will obtain its instructions from the transfer 
field as long as "hold transfer" (HT) is consecutively or- 
dered by another X punch in the order card. If the number 
occupying electronic storage C is positive at the time of the 
conditional transfer order, the computer continues to obtain 
its instructions from the normal instruction field. Alterna- 
tive routines can be programmed in the normal and transfer 
fields, and the discrimination between the two routines can 
always be phrased in such a way as to involve the sign of 
an appropriate quantity. 

If a given computation is to be performed for a series of 
different values of initially given parameters, it is necessary 
at each stage of the calculation to select the input parameters 
to be employed. It is frequently convenient to file with the 
order cards the input data necessary for a considerable num- 
ber of individual problems and to provide a means for the 
selection of that portion of the data required for a particular 



calculation. Provision is made for a selective loading opera- 
tion, termed "selective spread load" (SSL), controlled by 
an X punch in the order card. A load address, stored in the 
right-hand half of electronic storage D is compared with a 
load argument, punched in columns 14 to 18 of the order 
card. If the two numbers are equal, the X punch controls 
direct entry into the accounting machine counter groups 
2 to 7 in the same manner as in the ordinary spread load 
operation. If the two numbers are not equal, the card passes 
through the machine without activity. This order may also 
be employed in table lookup operations for the introduction 
into otherwise algebraic routines of transcendental func- 
tions, the direct calculation of which cannot be expeditiously 
programmed. The program controls are summarized in 
Table IV. 

TabIvE IV. Program Controls 



Punch 


Column 


Order 


Name 


X 


2 


CT 


conditional transfer 


X 


3 


HT 


hold transfer 


X 


5 or 35 


SL 


spread load 


X 


6 or 36 


SSL 


selective spread load 


9 


12 or 42 


P 


punch 


X 


12 or 42 


L 


list 



The CPC has been employed in this laboratory in the 
manner described on a number of different computational 
problems. The facilities for control have been found to be 
adequate for the efficient operation on computation routines 
of considerable complexity. We have also developed a set of 
control panels that are completely analogous, except that 
operations are limited to eight-digit arithmetic, and the 







45 146 47 48 49 50 51 52 53 54 55 50 57 58 59 (SO 61 52 63 64 «5 65 



Figure; 2. Controi, PANJii, for 417 Accounting Machine 
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square-root order is eliminated. In this way it is possible to 
preserve enough program capacity in the electronic calcu- 
lator for internal programming on the control panel of the 
604 unit for an iterative sub-routine for the calculating of 
transcendental functions such as the logarithm, exponential, 
trigonometric or hyperbolic trigonometric functions. 

We append notes on the various operations and list the 
sixty program steps for the general purpose ten-digit arith- 
metic control panel. A diagram (Figure 2) is given for the 
wiring of the accounting machine control panel, and it is 
supplemented with diagrams presenting, with greater clar- 
ity, calculate selector (Figure 3), field selector (Figure 4), 
co-selector (Figure 5), and pilot selector (Figure 6) cir- 
cuits. These notes are intended to supplement the material 
presented in the standard instruction manuals for the vari- 
ous machines. 



Notes 

Addition: On order 1, the operation A -\- B = C is per- 
formed. On order 6, the operation A -\- B -\- D = C is 
performed. 

Subtraction: On order 2, the operation A — B = C is 
performed. On order 7, the operation A — B -\- D = C 
is performed. 

Multiplication: On order 3, the operation ^ X ^ = C is 
performed. On order 8, the operation A X B -\- D = C 
is performed. The product is rounded off with 1/2 adjust- 
ment to ten significant figures. 

Division: On order 4, the operation A -^ B = C is per- 
formed. The formula C = (A—AiB2/Bi)/Bi is employed. 



-CALC.SEL- 

T O O 




(29) 



(28) 



(8) 
(28) 



(10) 



(20) 

(34) 

(19T^(32) 



(14) 



(30) 



(15) 



(1 

(2 

(3 

(4 

(5 

(6 

(7 

(8 

(9 

(10 

(11 

(12 

(13 

(14 

(15 

(16 

(17 

(18 



I.EGEND 

sup. without test (19 

sup. on zero (20 

PU group sup. 1 (21 

PU group sup. 2 (22 

PU group sup. 3 (23 
drop-out group sup. 4 (24 

program 3 (25 

program 10 (26 

program 13 (27 

program 15, 16 (28 

program 41 (29 

program 42 (30 

program 52 (31 

program 52 (32 

program 52 (33 

program 53 (34 

program 53, 55 (35 
read units into 2nd 



read units into 3rd 
read units into 4th 
read units into 5th 
read units into 6th 
read units out of 2nd 
read units out of 3rd 
read units out of 4th 
read units out of 5th 
read units out of 6th 
counter read in -|- 
counter read in — 
RI factor stor. 1 and 2 
RI general stor. 4 
RO factor stor. 1 and 2 
RO factor stor. 3 and 4 
RO general stor. 4 
RO mult. quot. 



SUPPRESSION CIRCUITS 



Group Suppression Uxits 



Other Suppression Circuits 



Suppression Program Lines Suppression 



Program Lines 



GS(1) 10-14 S (5) 6, 7,15,16 

GS(2) 4-5 S (6) 17,19,20,22,24,28,31,32 

GS(3) 8-9 S (7) 21 

GS(4) 18,23 S (8) 25,29,30 

S (9) 26,27,37-39 

S(10) 33-35,40,46-50 

S(ll) 36 

S(12) 41,42 

S(13) 43-45 

S(14) 52,53 

S(15) 54 

S(16) 55 



Figure 3. Calculate Selector Circuits of 604 Electronic Calculating Unit 




•(2) 

i< E G e; N D 
(1) Comparing Entry 1-5, Factor Storage 4 Entry (Split) (2) Comparing Entry 6-10, Multiplier-Quotient Entry (Split) 



FiGURt; 4. FiiCLD Selector Circuits oe 417 Accounting Machine 



.(1) 



(2) 



f COMPARING ENTRY / 

* m • — • • ' • ■ • •> • • 4 — » > 



o o o o o 



COMPARING EXIT 



ArVwG ENTRY >«-(4) ^^(3) 



o o o o o 




RBlIViiliriHU 



<HANNEL-A' 



■•rTTfnnn^lHlm 



lliliHilMiraL4J 



(1) field selector common 

(2) 2nd reading, cols. 14-18 

(3) pilot selector 8, IPU 

(4) general storage 4, exit 

(5) 2nd reading, cols. 34-38 

(6) 2nd reading, cols. 4-8 



(7) channel A control 

(8) field selector, PU 

(9) channel 13 control 

(10) 3rd reading, cols. 9-13 

(11) channel C control 

(12) shift control 



(21) 



LEGEND 

(13) digit selector 1, com. 

(14) operation control 

(15) 3rd reading, cols. 21-30 

(16) 3rd reading, cols. 31-43 

(17) 3rd reading, cols. 44-46 

(18) 3rd reading, cols. 47-50 



(19) 3rd reading, cols. 51-60 

(20) 3rd reading, cols. 61-62 

(21) 3rd reading, cols. 71-72 

(22) 3rd reading, cols. 14-20 

(23) 3rd reading, cols. 63-70 

(24) 3rd reading, cols. 73-80 

(25) pilot selectors 1-4, IPU 



Figure 5. Coselector Circuits oe 417 Accounting Machine 
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(1) split col. control 

(2) Co-selector 8, com. 

(3) card count 

(4) Co-selector 2, com. 

(5) comparing exit 

(6) field selector PU 



i.e;ge;nd 

(7) list 

(8) spread RI through col. split 

(9) spread RI and units position 
channel A control through 
col. split (split) 



(10) card cycles 

(11) set up change 1 

(12) 3rd reading, col. 2 

(13) 3rd reading, col. 3 

(14) 3rd reading, col. 42 

(15) 2nd reading, col. 12 



Figure; 6. Pilot Se;lEctor Circuit,-^ 
OF 417 Accounting Machine; 



where A^ = A (10-6), B^ = B (10-3), and B2 ^ B 
(2-1). The error is less than {A^/Bx) (^a/Si) 2. The left- 
most non-zero digit of B may not be more than one decade 
to the right of the leftmost non-zero digit of A. 

Square Root: On order 5, the operation \f'A = C is per- 
formed. The formula 

is employed, where C^"^ is the «th approximation to the 
square root. The successive approximations are continued 
until for some n, say n', C^"'^ = C^^-'+^K The zeroth ap- 
proximation is obtained from the empirical rule, 

C(«> = 2(10^-i«+«/2) -f yi( 10^-1 +«/2) , 



where the number of non-zero digit pairs of the 10-digit 
number ^ is ^ -|- 4, and where a denotes the location of the 
decimal point of A, counting from the leftmost decade of A. 
The quantity A must be transmitted over channel A at the 
time square root is ordered. It is also essential that a be an 
even number. 

Transfers: The C-^A transfer accompanies every arith- 
metic order. The C-^B and C—>D transfers can be ordered 
with every arithmetic order except square root. The D-^B 
transfer can accompany simple addition, subtraction, multi- 
plication, and division (orders 1, 2, 3, 4) only. 

Spread Load: Counters 2 to 7 read in from card fields. 
Counter 1 can be loaded simultaneously by card entry into 
channel A if channel C is given a 71 address. 
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Selective Spread Load: Counters 2 to 7 read in from 
card fields of the card for which the argument (cols. 14-18) 
is equal to a card address stored in electronic storage 
i>2 = D (5 — 1 ) . All cards for which the argument and card 
address are non-equal are passed through the machine with- 
out reading to the counters. 

Conditional Transfer: The calculator transfers to the 
transfer instruction field instead of the normal instruction 
field if C < when conditional transfer is ordered. Hold 
transfer must also be ordered. 



Hold Transfer: The calculator reads instructions from 
the transfer instruction field instead of the normal instruc- 
tion field if transfer has occurred as a result of conditional 
transfer, as long as hold transfer is continuously ordered. 
The first transfer instruction word must be blank, and the 
first normal instruction word after transfer must be blank. 
If either normal or transfer instructions call for entry into 
channel A or B from card fields, the opposing instruction 
word must be blank. 



TYPE 604 PROGRAM STEPS FOR TEN-DIGIT ARITHMETIC 





Sup 


Sel 


1 
Factor Storage 


Mult 
Quot 


Ctr 


General Storage 


Activity 


No. 


1(8-6) 
2 


3{8-6) 
4 


1(8-6) 
2 


3 


4 


+ 


— 


X 


•4- 


r 


1 







(PUGS 


4) 




RO,RC 








X 


X 


X 


X 


X 


2 













RI(+)2 




RO 




X 


X 


X 


X 


X 


3 













RO,RC 




RI 




X 


X 


X 


X 


X 


4 


GS2 




R14 


RO 


















X 


X 


5 


GS2 




(PUGS 


4) 


RO 








RI 








X 


X 


6 


5 






RO 




RI(+)4 








X 


X 






X 


7 


5 






RI 




RO,RC 








X 


X 






X 


8 


GS3 






RO 




RI(+)6 








X 


X 






X 


9 


GS3 








RO 


RI(+)4 








X 


X 






X 


10 


GSl 


IN 




RI 




R06,RC 
















X 






IT 




RI 




R05,RC 
















X 






2T 




RI 




R04,RC 
















X 






3T 




RI 




R03,RC 
















X 






4T 




RI 




R02,RC 
















X 


11 


GSl 






RO 




RI(+)2 
















X 


12 


GSl 




(EMITTER 4) 


RI6 


















X 


13 


GSl 


IN 






RO 


RI(+)5 
















X 






IT 






RO 


RI(+)4 
















X 






2T 






RO 


RI(+)3 
















X 






3T 






RO 


RI(+)2 
















X 






4T 






RO 


RI(+) 
















X 
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Sup 


Set 


Factor Storage 


Mult 
Quot 


Ctr 


General Storage 


Activity 


No. 


3(8-6) 
2 


3(8-6) 
2 


1(8.Q) 
2 


3 


4 


+ 


_ 


X 


-1- 


V 


14 


GSl 










R03,RC 


RI 














X 


15 


5 


2N 
5N 








RI(+)6 


RO 






X 














2T 

5N 








RI(-)6 


RO 








X 












5T 








RI(+)6 


RO 














X 


16 


5 


2N 
5N 








RI(+)3 




RO 




X 














2T 

5N 








RI(-)3 




RO 






X 












5T 








RI(+)3 




RO 












X 


17 


6 






RO 




RI(-)5 
















X 


18 


GS4 








RO 


RI(-)3 
















X 


19 


6 










R03 




RI 












X 


20 


6 










R06 


RI 














X 


21 


7 




(PROGRAM R] 
(ZERO. TEST) 


3PEAT) 


R16 
















X 


22 


6 










RO.RC 
















X 


23 


GS4 








RO 


RI(+)2 
















X 


24 


6 




1 
(EMITTER 5) 


RI4 


















X 


25 


8 










MULT(+) 


RO 










X 




X 


26 


9 






RO 


RI 














X 






27 


9 






RI 




R06,RC 












X 






28 


6 




(EMITTER 5) 


RI 


















X 


29 


8 










MULT(-I-) 




RO 








X 




X 


30 


8 






RO 




RI(+)4 












X 




X 


31 


6 








. 


R02 




RI 












X 


32 


6 










R05,RC 


RI 














X 


33 


10 




(PUGS 


1) 


RI 






RO 










X 


X 


34 


10 




RO 






MULT(+) 














X 


X 


35 


10 




(PUGS 3) 


DIV 




RO 












X 


X 


36 


11 






RI 




R04,RC 












X 


X 


X 


37 


9 










MULT(-I-) 


RO 










X 
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Type 604 Program Steps for Ten-Digit Ari 


THMETIC- 


-Continued 






Sup 


Sel 


Factor Storage 


Mult 
Quot 


Or 


General Storage 


Activity 


No. 


3{8-6) 
2 


3{8.6) 
2 


1(8-6) 
2 


3 


4 


+ 


— 


X 


-;. 


V 


38 


9 






RO 




RI(+) 












X 






39 


9 




(H ADJUST) 




RI3 












X 






40 


10 








RO 


RI(-)2 














X 


X 


41 


12 


6N 


RO 






RI(+)6 














X 


X 






6T 


RO 






RI(+)6 








C+D 








42 


12 


6N 








RI(+)3 






RO 








X 


X 






6T 








RI(+)4 






RO 


C+D 








43 


13 








RI 


R04 








X 


X 


X 






44 


13 








RO 


RI(-)4 








X 


X 


X 






45 


13 






RI 




R06,RC 








X 


X 


X 






46 


10 








DIV 




RO 












X 


X 


47 


10 






RI 




RO.RC 


(PUGS 2) 










X 


X 


48 


10 






RO 




RI(+)6 














X 


X 


49 


10 






RI4 
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Remarks on Distillation Calculations'^ 

JOHN W. DONNELL 

Michigan State College 



M 



THE WORK we have been doing with IBM calculating 
machines has been in distillation and absorption tower 
design. 

First, I should like to say a few words regarding Mr. 
Opler's paper^ on distillation and the previous paper by 
Dr. Rose^ and his associates, to avoid repetition in describ- 
ing our work which covers similar ground to that of these 
two excellent papers. It should be pointed out that where 
Mr, Opler assumed Raoult's law to hold, he is using the 
equation 3; = PX/ir to describe the relationship between the 
concentration of a given component in a boiling liquid and 
the concentration of that component in the vapor in equi- 
librium with the boiling liquid. In this equation, tt is the 
total pressure on the system, P is the vapor pressure of the 
component in question in the pure state at the temperature 
in question, and x and y are the mol fractions, respectively, 
of the component in the liquid and vapor states. 

In our work, especially in petroleum engineering, we 
have discovered that such an equation is not accurate 
enough for our design work, and the petroleum industry 
has quite universally adapted the equation 3; = KX where, 
as you can see, P/tt has been replaced by K, an experi- 
mentally determined factor, which is a function of tempera- 
ture, pressure, and the component itself. These constants 
have been determined at great expense to the industry, and 
I think the industry should be called upon to pool its K 
factors. I suggest that IBM make a master deck of cards 
containing these values. In our work we are using K factors 
in both distillation and absorption calculations. In Dr. 
Rose's paper a symbol designated as a is, in reality, a ratio 
of K values; e.g., ai_2 = Kx/K^. While this varies much 
less with temperature than the individual K values, we in 
the petroleum industry cannot use this approximation in 
most design work, and must use the individual K factors 
and make no attempt to eliminate variations with tempera- 
ture in either distillation or absorption columns. 

*This paper was presented informally. 



In addition I should like to point out the similarity, in 
the approach to the distillation problem, o£ the method we 
are using to that of Dr. Rose. You recall his calculations 
concern themselves only with binary mixtures. He calcu- 
lates the analysis of the distillate and bottoms products for 
a given feed, reflux ratio, and number of plates, by first esti- 
mating the total amount of distillate per unit feed and then 
calculating a new distillate quantity closer to the correct 
answer. With this new distillate quantity as a second esti- 
mate, he calculates a still closer answer and continues the 
iteration until the desired accuracy is obtained. We do a 
similar calculation except that we do not assume a constant 
a and we work with multicomponent mixtures. After ob- 
taining our value of the amount of distillate by iteration, we 
must go back and correct for variations in K values caused 
by our inaccurate estimation of the correct temperature at 
the start of the first iteration. This, then, imposes a second 
iteration upon the first. Although this may sound compli- 
cated, it is really not so difficult, because we have prepared 
a detailed calculation chart showing just how each calcula- 
tion comes from another (Figure 1). This chart is simple 
enough that one not skilled in the art — say a high-school 
student — could carry the calculations out on a desk calcu- 
lator. Also, we have prepared a chart for absorption which 
we have adapted to the IBM Type 602 Calculating Punch, 
and we are in the process of doing the same for the distilla- 
tion chart. In using the type 602 machine we have to use 
some seven control panels (wired separately) to carry out 
the various calculations. We are also considering, with the 
help of the IBM Applied Science Department, a somewhat 
simple and direct attack on the distillation calculation, using 
matrix algebra to solve a large number of equations used to 
describe the conditions in a distillation tower. 

re:fe:ri:ncss 

1. AschEr OplEr, "Machine Calculation of the Plate-by-Plate 
Composition of a Multicomponent Distillation Column," pp. 18-23. 

2. Arthur Rose, Theodore J. Williams and William S. DyE III, 
"Continuous Distillation Design Calculations with the Card-Pro- 

„ grammed Electronic Calculator," pp. 24-31. 
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Figure 2. Absorption Cai.cui.ations 



Some Applications of the Nlonte Carlo Nlethod 



NLatrix Inversion on the 
IBM Accounting Machine^ 

ASCHER OPLER 

The Dow Chemical Company 

IN THE JULY, 1950, issue of Mathematical Tables 
and Other Aids to Computation, Forsythe and Liebler de- 
scribed a Monte Carlo method of matrix inversion. This 
has been adapted for use with an IBM Type 405 Account- 
ing Machine. In the interest of brevity, I state a version of 
the original method without proof. 

If 5 is a square matrix such that A = \ I — B \, then 



B- 



CO 

{\ I — A\)-'^ = \ A^. When the last expression 



is a convergent series, B~^ may be evaluated by playing an 
infinite number of games as described below. In practice it 
may be approximated by playing a large number. 



Choose stop probabilities, pj, such that 



w 



+ PJ = 1- 



(This is one restriction of the class of matrices that may be 
inverted.) Prepare a deck of punched cards as follows: 
one half of the card is to contain random numbers; a playing 
field of n columns is prepared in the other half of the card. 
In each of the n columns, punch digits so that the proba- 
bility of a card containing a^ or pj is equal to the correspond- 
ing element of A or the corresponding stop probability. 
"Shuffle" by repeated sorting in the randomizing field. 

The game is played by passing this prepared deck 
through the accounting machine n times. To play for in- 
verse element (Bij)~^, the machine is instructed to start 
with the first card and read the ith. column on that card. If 
the number read is m, the machine reads the wth column of 
the next card. Each card directs the reading of the following 
card. When a stop card is read, the game is ended. If the 
stop card was read in column ;', the score is 1; otherwise the 
score is zero. Actually, n games are played at once with one 
of the n elements winning and others losing. The inverse 
elements are found by dividing the final scores by total 
games played and then finally dividing by the stop proba- 
bility. In certain cases, the inverse matrix may be printed 
by the accounting machine, line-by-line. 

*This paper was presented informally. 



This method appears to be the simplest and most econom- 
ical way of inverting matrices. The operations are pro- 
portional to n^ instead of n^. The disadvantage is that the 
results obtained in reasonable times are approximations. 
(A fair approximation to a seventh-order inverse may be 
found in less than an hour; a good approximation might 
take four to six hours.) The class of matrices that, may be 
inverted is limited, but, with ingenuity, many of these limi- 
tations may be overcome. 

A paper describing the method fully has been submitted 
for publication in Mathematical Tables and Other Aids to 
Computation. 



Remarks on Finding Roots ofj 
and Inverting^ a Matrix^ 

GILBERT W. KING 

Arthur D. Little, Incorporated 

Ths ESSENTiai, feature of the Monte-Carlo Method of in- 
verting a matrix, or of finding its roots, is the well-known 
iterative procedure of raising the matrix to a high power. 
The reason the method does this can be easily seen as 
follows : Choose a row and column at random and write 
down the matrix element, ai\. Then choose another element 
from the Ath row lying in column /x, say, axu,. Again choose 
an element from the /nth row, say a^y. This is done N times. 
The choices are multiplied together, 

fliX a\n a^v • • • O'Tij • 
We recognize this as a term in the ijth. element of the Nth 
power of the matrix. If we took all possible choices of paths 
from row to column, starting at the ith row and ending a 
yth column, and added the products, we would have pre- 
cisely the value of ijth. element of the Nth power. The 
random procedure described above merely picks some terms 
from some elements at random. By having the probability 
of picking any element, say, a^^, proportional to its magni- 
tude, the Monte-Carlo Method picks out the terms in pro- 
portion to their magnitude, and hence gets the principal 
term of the principal elements of the Nth power of the 
matrix.^ 

REFERENCE 

1. See also Gilbert W. King, "Stochastic Methods in Quantum 
Mechanics," Seminar on Scientific Computation, November, 1949 
(IBM). 

*This paper was presented informally. 
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Remarks on the NLonte Carlo Method 

CUTHBERT C. HURD 

International Business Machines Corporation 

The: Month Carlo Method has aroused great interest 
and has found many applications. A general description of 
the method as applied to a class of problems in mathemati- 
cal physics is given by Metropolis and Ulam.^ The Proceed- 
ings of a Symposium on Monte Carlo Methods, which was 
held in Los Angeles from June 29 to July 1, 1949, is to be 
published by the National Bureau of Standards. Finally, 
papers on the subject have been presented at the IBM Sem- 
inar on Scientific Computation, November, 1949,^ and the 
IBM Computation Seminar, December, 1949.^ 

In applying the method it is convenient to have available 
a set of random digits. These digits can be computed by 
several schemes which have been proposed. For example, 
60,000 such digits can be computed in an hour on the IBM 
Type 604 Electronic Calculating Punch. Once the set of 
random digits is available, sets with a prescribed probability 



distribution can be obtained by performing a table lookup 
operation using the sorter and collator, or in some cases by 
direct calculation. 

The Monte Carlo Method is useful in many instances in 
giving an estimate of an answer. The results of mathemati- 
cal statistics then allow one to attach a measure of reliability 
to the answer. 

REl^eRENCES 

1. N. Metropolis and S. Ulam, "The Monte Carlo Method," Jour. 
Am. Stat. Assoc, Vol. 44 (1949), pp. 335-341. 

2. Seminar on Scientific Computation, November, 1949 (IBM). 
William W. Woodbury, "Monte Carlo Calculations," pp. 17-19. 
Herman Kahn, "Modification of the Monte Carlo Method," 
pp. 20-27. 

Z. W. BiRNBAUM, "On the Distribution of Kolmogorov's Statis- 
tic for Finite Sample Size," pp. 33-36. 

Gilbert W. King, "Stochastic Methods in Quantum Mechanics," 
pp. 42-48. 

John H. Curtiss, "Sampling Methods Applied to Differential 
and Difference Equations," pp. 87-109. 

3. Computation Seminar, December, 1949 (IBM). 

Mark Kac and M. D. DonskEr, "The Monte Carlo Method and 
Its Applications," pp. 74-81. 

P. p. Johnson and F. C. Ueeelman, "A Punched Card Appli- 
cation of the Monte Carlo Method," pp. 82-88. 
Everett C. YowEll, "A Monte Carlo Method of Solving La- 
place's Equation," pp. 89-91. 

Gilbert W. King, "Further Remarks on Stochastic Methods in 
Quantum Mechanics," pp. 92-94. 



Plotting Punched Card Data Using the 
IBM Type 405 Accounting Machine"^ 



PAUL T. NIMS 
Chrysler Corporation 



A METHOD for plotting the points on a curve has been 
described in the writings^ of Dr. Gilbert W. King. This 
method gives a plot having abscissa values from 1 to 80 and 
an unlimited range of ordinate. However, for the purpose 
of checking some automotive design calculations, it is con- 
venient to have a larger range of ordinates and a relatively 
small range of abscissas, together with adaptability to mul- 
tiple-valued functions. 

A method similar to that of Dr. King's but differing in 
some particulars was worked out by Mrs. Virginia Johnson 
and Mr. F. F. Timpner of the Chrysler Corporation Engi- 
neering Division. First, a master deck was prepared, cover- 
ing the range of ordinate values desired (positive and nega- 
tive) with the ordinates punched in columns 79-80 as well 
as in the same field as the ordinates on the problem cards. 
The master deck was distinguished with an X punch in 
some suitable column. 

The master deck and the problem deck were sorted to- 
gether (in that order) on ordinate values to give a com- 
bined deck ready for the accounting machine. 

As mentioned above, the range of abscissas is relatively 
small; in fact, determined by the available selector and type 
bar capacity, the range is to 69. The two columns are read 
by the 405, and the tens and units digits enter separate digit 
selectors. The tens digits (one through six) control the six 
different 10-position field selectors on the type 405 machine 
to pick up one of them, depending on the tens digit read. 

The ten units digits are wired, position for position, to 
the 10 common hubs on the field selector 6 and from its 
transferred hubs to counter input hubs 60-69. The normal 
hubs of field selector 6 are wired to the common hubs of the 
field selector 5, the transferred hubs of the field selector 5 
are wired to counter inputs 50-59, and the remaining selec- 
tors are laced similarly. The normal hubs of the field se- 
lector 1 are wired to the first ten counter positions to take 
care of abscissas through 9. The only exception is abscissa 

*This paper was presented informally. 



43, which is not wired to any counter, as the 405 has a type 
bar support at this point. The counter exits are wired, posi- 
tion for position, to the total type bars. 

Counters are all wired to add from plug to C, but carry 
wiring is unnecessary. The accounting machine is wired to 
print a minor total on the X in the master cards. Alterna- 
tively, if the ordinate values are punched in the same col- 
umns on both detail and X cards, a minor total pickup can 
be used. This minor total, also, advances the paper one line. 
With this arrangement, the abscissas of all cards having the 
same ordinate are stored in the counters, and the whole 
line is printed at once. All of the zero suppressors are raised. 

There are a few special wires. Columns 79-80 of the 
master cards are listed to give an ordinate scale. As no 
zeros will print, the abscissa values ending in zero are wired 
to an X selector, and a card count pulse is entered in the 
counter in order to print a 1 for the zero values. A set of 70 
list wires is also used to print headings on the curve and to 
print the abscissa scale. 

This plotting method is very convenient for checking the 
results of engineering computations, and the illustration 
(Figure 1) shows a function which had a few deliberate 
errors introduced in the last place to check the sensitivity of 
error detection. 

reb^icre^nce: 

1. Gilbert W. King, A Method of Plotting on Standard IBM 
Equipment, MTAC, Vol. 3, No. 25 (January, 1949), pp. 352-355. 

DISCUSSION 

Mr. Bejarano: What is the possibility of using the reverse 
process and picking up the number from the graphs? 

Dr. Hurd: Various systems have been tried and pro- 
posed. I will tell you what I know about it. 

Mr. Bejarano: I have in mind distillation, using the 
graph of the case where the charts agree without any arti- 
ficial function. 
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Dr. Hurd: One device, which is available at the present 
time, which goes from a plot to IBM cards is that which is 
done by manual operation to locate the ordinate and the 
abscissa. In fact, there are at least two devices of that kind. 
There has been experimentation on using photoelectric 
cells. 

Dr. King: Some of the charts on the slides I have were 
made on the IBM accounting machine. We did not con- 
sider it a trivial proposition. We had hundreds of these 
plots to do. In our scheme we don't have any sorting and 
collating to do. We take the cards as they come out of the 
computer; the ordinate this way can be in any field of the 
card. Just change a few of the wires, and, as the cards feed 
through the machine, points are plotted just as fast as list- 
ing speed. There are no total cycles. 

Dr. Hurd: Many variations of these schemes have been 
worked out. One scheme used at the Boeing Aircraft plant 
was to make a mark which is an 8. They also put in calibra- 
tion marks and the expected average mark. Then they 
photograph it and reduce it, by which time these 8's look 
like dots and the graph looks like a continuous graph. You 
have essentially two-digit accuracy of plots here. As the 
cards are run through, the first differences are calculated 



and printed, so that if you want to interpolate, you have at 
least first order interpolating coefficients already plotted on 
the sheet. 

Mr. Mongreiff: All of these schemes depend on selecting 
correct type bars. This is especially useful where the length 
of the graph is much greater than the height, but if you 
could turn the system around, that is, if your needs are that 
the height of the graph is greater than the length, it seems 
to me you could depend on successive reduction on the total 
in a counter, put the data in a counter to begin with, and 
then reduce it by one or two or three successively. You 
could even make that amount manually changeable just by 
shifting a wire or two and then emitting an X or another 
impulse when the counter turned negative. It would be 
possible to start off with a complement, say a 100 comple- 
ment or a 500 complement of the original. In that way you 
would print a graph which was right side up. 

Dr. Hurd: The 407 is a new accounting machine which 
we have not demonstrated here. Perhaps that would be 
capable of doing this. 

Mr. Opler: In mentioning the 407, it reminds me that the 
407 would be twice as powerful a tool for graphing as the 
402, because it has 120 type bars across and lists 150 lines 
per minute. You can interpret in a hurry. 



A JS/lethod for Evaluating Determinants and 

Inverting Nlatrices with Arbitrary Polynomial Elements 

by IBM. Punched Card Methods 



L. E. GROSH, JR. 

Purdue University 
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THE METHOD of determinant evaluation proposed 
here is essentially a scheme for the collection of terms with 
like coefficients and like powers of x, where these terms are 
polynomials with non-numerical coefficients. The method is 
based upon the termwise expansion of a determinant; thus, 
it is concerned with the evaluation of n ! polynomials for an 
nth order determinant. If any zero elements appear in the 
determinant, a significant reduction in the number of cards 
handled may be made. The method works best when there 
are few distinct coefficients in the original determinant. 

The inversion of a matrix with polynomial elements is 
based upon the expression for the general element of the 
inverse which involves the ratio of two determinants. Thus, 
for an wth order matrix, the inversion problem is reduced 
to the evaluation of «^ determinants of order (m— 1) and 
one wth order determinant. Most of the cards necessary for 
the evaluation of the n^ determinants are generated in the 
evaluation of the wth order determinant. 

The equipment needed for the evaluation procedure is: 
a card punch, sorter, reproducer, 602-A calculating punch, 
and 405 accounting machine — direct subtracting. 
The use of the calculating punch is not absolutely necessary; 
its role in the procedure could be performed by a sorting 
and gang punching operation. However, for a problem of 
any size, the sorting operation would become extremely 
complicated and impractical unless an IBM Type 101 Elec- 
tronic Statistical Machine were available. 

EVAI^UATION 01'' DETE^RMINANTS 

The Coding Problem 

The heart of the evaluation scheme is a method of coding. 
The coding is best understood when applied to a simple 
example. Consider the third order determinant shown below. 



The termwise expansion of this determinant is: 

D = Pll{x)P22{x)Pzz{x) + P21{X)PUX)P^^{X) 

+ Pzi{x)Px2{x)P2^{x) -P3l(x)P22ix)P^s{x) 
- P2lix)Pr2{x)P^six) - Pll{x)Ps2(x)P2zix) 

If the factors of each term are ordered on the column sub- 
script as above, these six terms may be represented by the 
permutations of the numbers 1, 2, 3, and an X punch to 
indicate a negative sign. Thus, 

123 = Piiix)P22{x)Ps3(x) 

231 =P2l(^)i'32(^)-Pl3(^) 
m = P3^(x)P,2ix)P2B(x) 

321X = -Pzx{x)P22{x)Pi3{x) 
213X = -P2y{x)Px2{x)Pz3{x) 
U2X= -Pir(x)Ps2(x)P2s(x) . 

Let each one of the six symbols 123, 32 IX, ... be called a 
permutation number. 

Now consider one of the products, say 
P2x(x)Pi2(x)P33(x) = 

(al^x + a^i) {al^x + aj^) (afx + af ) 

= af^ a\^ af^x^ 

+ « of af + af ai2 af + a^i al^ af )^2 

+ al' af af . 
Note that there are 2(2) (2) = 8 different terms in the 
expanded polynomial. In general, there will be at most, 

n 

TT {kij +1) terms in this expansion, and there will be 

exactly this number if all o^ =7^ and are distinct, where kij 
is the degree of Pij{x), and n is the order of the determi- 
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nant. Each of these three-factor coefficients may be repre- 
sented by a three-digit number — namely, the subscripts 
of the a%'s involved. Thus, 

111 = afi af af 100 = af aj^ af 



110 = ofi al^ af 
101 = a?i a 



Oil = agi a 






010 = a^i al^ af 
001 = a2i aJ2 af3 

000 = a21 ^12 ^33 



Let each of these three-digit numbers be called a selection 
number. The first digit i of the selection number ijk refers 
to the coefficient of .«■* from a polynomial in the first column; 
the second digit / refers to the coefficient of x' from a poly- 
nomial in the second column; and the third digit k refers to 
the coefficient of x^ from a polynomial in the third column. 
Also, the sum of the digits i-\-j-\-k is the power of x with 
which this three-factor product is associated. 

The combination of a permutation number and a selection 
number allows us to code each term in the expansion of D; 
also to indicate its sign and the power of x in the expansion 
with which the three factor coefficient is associated, Thus, 



if 



biX^ -\- &2^ + &3 



132 Z 302= -al^ a^ » 
and is a coefficient of x^ in the final expansion. 

If all of the a%'s are distinct, the coding problem is fin- 
ished, but in some important cases not all of the a^'s will be 
distinct; thus there is one more coding step to be consid- 
ered. For example, take the following special case of the 
above third order determinant D. 

bsx^ -f- b2X^ + ^1-^ + ^0 b^x -\- &2 
bs^ + bo b2X -|- bi 

&2 bix -f bo 

There are many ways in which the coefficient ± &o^2&3 may 
arise, e.g., 

132X200= -&2&0&3 

132 Z 301 = -bsbobz 

213 Z 101 = -&3M0 

213 Z 000= -bobzbs 
etc. 

There are only four distinct coefficients in our problem; 
therefore a four-digit number will be sufficient to identify 
any possible combination of these biS. These will be called 
term numbers. Let the first digit of this number be the 
power of bo in the term, the second digit be the power of &i, 
etc. Thus, 

1011 = b^bobs 
= &3&0&2 
= &3&2&0 

= &0^2&3 . 

The main points of the coding scheme can be summarized 
as follows: 

1. A three-digit permutation number indicates which 
three polynomials are multiplied together and indi- 
cates the sign of the resultant polynomial. 



2. A three-digit selection number indicates which co- 
efficient of each polynomial shall be multiplied to- 
gether to obtain a resultant coefficient. The power of 
X, with which this resultant coefficient is associated, 
is the sum of the digits. 

3. A four-digit number indicates the power of the origi- 
nal distinct coefficients as they appear in the resultant 
coefficient. 

Card Generation 

Although this coding scheme is relatively simple, it 
would be almost impossible to maintain 100% accuracy on 
even a medium-size problem if the coding were done by 
hand. Fortunately, all of the coding may be done by ma- 
chines with very little key punching. To accomplish this, 
four different decks of cards are used. These are: 

1. Permutation deck 

2. Selection deck 

3. Master gang punching deck 

4. Term deck 

The term deck is the final product, and the other three decks 
are used in its generation. A typical card layout for the 
problem considered above is as follows: 



Card Layout 



Permutation 
Deck 



Selection 
Deck 



col. Term Deck 

2 ^ permutation number 

3J . .... 

4 X if permutation number is positive 
X if permutation number is negative 



5 

6 

7 

8 

9" 
10 
11 
I2I 
13 
14 
15 
16 



selection number 

the actual coefficients indicated by the 
permutation and term numbers 

bo=0 &i = l &2 = 2 &3 = 3 

term number representing the power of 
the various bm's as they appear in the 
' final answer. 



the power of x with which the term 
number is associated (the sum of col- 
umns 6, 7, and 8). 
17 a check of the machine operation, the 
sum of columns 12, 13, 14, and 15. For 
this problem it should be a 3 for every 
term card. 
The master gang punching deck will have a distinguish- 
ing X punch, say in column 80, and is divided into three 
parts: 

1. punching in columns 1, 6, and 9 

2. punching in columns 2, 7, and 10 

3. punching in columns 3, 8, and 11. 
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This deck is used to facilitate the punching of columns 9, 
10, and 11 of the term deck after the first eight columns 
have been generated from the permutation and selection 
decks. Consider part 2 of this deck; the required cards are: 

Column 
2 7 10 



1 3 = &3 

1 1 2 = &2 

2 2 = &2 

2 1 1 = &i 

3 1 = bx 
3 I = bo. 

The generation of the permutation and selection deck is 
quite simple, but is a very important part of the coding oper- 
ation. The permutation deck for an nth. order determinant 
is generated from the permutation deck for a determinant of 
the («— l)st order. If the generation scheme is carried out 
correctly, it is only necessary to handle n ! cards for this 
operation. While it is possible to generate the permutation 
deck for any order determinant from a single key punched 
card, it is more practical to select some moderate-size order, 
say 4th, key punch the required 24 permutation cards, and 
start the generation procedure from this point. Since the 
example considered above is of 3rd order, we shall start 
with a 2nd order determinant for illustration. For ease of 
machine operation we shall use two columns to designate 
the sign of the permutation. These columns shall be the 
(w-)-l)st and the («-)-2)nd columns, if we are generating 
cards for an nth. order determinant. Thus, for our example, 
the sign will be in the 4th and 5th columns. If the generation 
procedure is started with permutation cards for an mth 
order determinant and cards for an «th order determinant 
are desired, the following steps must be repeated (n — m) 
times: 

1. Letting — denote a blank column, the permutations 
for a 2nd order determinant are 

12-X~ and 21-- -X, 
the X in column 4 denoting a positive permutation 
and the X in column 5 denoting a negative permuta- 
tion. 

2. Gang punch a 3 in column 3, giving 

123X- and 213 -X 

3. Reproduce the cards from step 2, interchanging col- 
umns 2 and 3 and columns 4 and 5. Thus, 

123X- gives 132 -X, 
213-X gives 231X-. 

4. Reproduce the cards from step 3, interchanging col- 
umns 1 and 2 and columns 4 and 5. Thus, 

132-X gives 312X-, 
231X- gives 321 -X. 

5. For larger n, this procedure of interchanging columns 
is repeated until the n gang punched in step 2 in the 
Mth columns has been moved over to the first column. 



Thus, there are («— 1) interchanges. 
The six cards we have generated are: 

123X- 132-X 312X-- " 
213-X 231X- 321-X, 
and are all the permutation cards necessary for a 3rd order 
determinant. 

The steps in the generation of the selection cards are: 

1. Determine the highest degree appearing in each col- 
umn of the determinant. The maximum number of 

n 

cards needed will heTT (kij-{-l) as stated above, 

where kij is the highest degree appearing in column /. 

2. Using blank cards, key punch the numbers 0, 1, ... , 
kii into column «-f 3, one number to each card. For 
our example this would be 0, 1, 2, and 3 in column 6. 

3. Gang punch a in column «-|-4 of the cards made in 
step 2, giving 00, 10, 20, and 30. 

4. Reproduce column w-|-3 of step 3 and gang punch 1 
into column n-\-4, giving 01, 11, 21, and 31. Repeat 
this operation until ^2/ has been gang punched into 
column w-f-4. 

5. Gang punch a in column n-\-5 of the cards gen- 
erated in steps 3 and 4. 

6. Reproduce columns (n-\-3) and (n-\-4) and gang 
punch a 1 in column («-f-5). 

7. Repeat step 6 until k^j has been gang punched into 
column (w-f5). 

8. The operations of gang punching and reproducing are 
repeated until knj has been gang punched into column 
(2«+2). 

The 24 cards generated for the determinant D by the above 
procedure are: 

000 010 100 110 200 210 300 310 

001 Oil 101 111 201 211 301 311 

002 012 102 112 202 212 302 312 

The next step in the evaluation procedure is the genera- 
tion of the term cards. The general procedure to be followed 
here is to place the selection cards in the reproducing hopper 
and a blank deck of cards preceded by a single permutation 
card in the gang punch hopper. Then, gang punch the per- 
mutation number and reproduce the selection number into 
the blank deck, column-for-column. This procedure would 

n 

yield n ! j 1 (kij-\- 1 ) term cards, many of which would have 

no meaning unless all elements of the same column had the 
same degree. Referring to the example above, 3 !(24) == 144 
term cards would be generated by the unmodified repro- 
ducing and gang punching procedure. By a study of the 
degree of the polynomial represented by each permutation 
card, it is easily seen that only 60 term cards are needed 
for the above example. 
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Permutation 




Number 


kxi^ 


123X- 


4 


132 -X 


4 


213-X 


2 


231X- 


2 


312X- 


1 


321 -X 


1 



Number of 
k2j-\- 1 ksj-j- 1 Term Cards 



16 

24 

8 

4 

6 

2 



Total number of term cards 60 

The extraneous term cards can be eliminated and much 
time and many cards saved if a sorting procedure is used 
before the term cards are punched for each permutation 
card. As an aid to the sorting operation, rewrite the de- 
terminant in the following symbolic form: 

0, 1, 2, 3 0, 1 

0, 1 0, 1 0, 1, 2 

0, 1 0, 1 

where each element has been replaced by the subscripts of 
the o^ s appearing in that element. From this symbolic form 
it is easily seen which selection numbers are needed for any 
permutation number. The sorting procedure to be followed 
for each permutation number is: 

1 2 



Permutation 
Number 



Sort on col. 6 
Reject 



Sort cards from i 

on column 8 

Reject 



123X- 
132 -X 
213-X 
231X- 
312X- 
321 -X 



none 
none 
2,3 
2,3 
1,2,3 
1,2,3 



2 
none 

2 

1,2 
none 
1,2. 



The selection cards left after this sorting procedure are the 
ones to be used in punching the term cards. 

Machine Operations 

After the generation of the term cards, the following ma- 
chine operations must be performed: 

1. Using the master gang punch deck, gang punch in the 
appropriate columns (9 to 11) the coded values of 
the bi's. 

2. Using the 602-A calculating punch, count the number 
of times the various biS appear in columns 9 to 11 
and punch this information in columns 12 to 15. 

3. On the 602-A, crossfoot columns 6, 7, and 8 to deter- 
mine the power of x with which each term card is 
associated; punch this sum in column 16. At the same 
time, crossfoot columns 12, 13, 14, and 15 and punch 
in column 17 to obtain a check on operation 2. 

4. By sorting, arrange the cards in order by term num- 
ber, columns 12 to 15, within each power of x, col- 
umn 16. 



5. Using a 405 or any other direct subtraction account- 
ing machine, tabulate the term cards, controlling on 
the power of x and the term number, columns 12 to 
16, adding a 1 for each positive term card and sub- 
tracting 1 for each negative term card. Either column 
4 or 5 may be used for this purpose. Naturally, the 
term number and power of x should be group indi- 
cated. 

Inversion of Matrices 
Using the expression 



A, 



an 



det(aiy) 

for an element of the inverse of the matrix {aij), where Aij 
is the co-factor of a^, it is easy to adapt the above determi- 
nant methods to evaluate this expression. The evaluation of 
det(aiy) needs no explanation; only the cofactor A^ must be 
considered here. We can write 



det (Oi;) 



— 7 aij/iij . 



To evaluate A^, it is only necessary to select from the term 
cards used in the evaluation of det(ai^-) those cards which 
contain the element an. If these cards are tabulated in the 
same manner as for det(aij), the result will be a^Aij, and to 
obtain A^ we must factor out a^- from this tabulation. Be- 
cause we are considering polynomial elements, 

aij = rij\X) , 

the factoring procedure would be cumbersome. To avoid 
this, instead of selecting all of the term cards for Pij{x), 
select only those for a particular a'^ of Pi}{x). Then a^can 
be factored out of the tabulation record by subtracting 1 
from each entry in the column of the record correspond- 
ing to a^ which can be done automatically by the accounting 
machine itself. 

This technique will not give us all the cofactors we need 
if some of the a^'s in the original matrix are zero, because 
in the evaluation of det (aiy) we have not punched the term 
cards corresponding to the a^y = 0. In the case where there 
are some a^ = 0, the following procedure may be used in 
the generation of term cards: 

1. Select all ay ^ permutation cards as before. 

2. Of the permutation cards with at least one ay = 0, 
select those in which exactly one ay = appears. 

3. Generate the term cards as before, except that in the 
master gang punching step use a 12 punch to identify 
those term cards corresponding to an ay = 0. 

4. Use these new cards to obtain Aij corresponding to 
ay = 0. It is no longer necessary to factor out the ay 
because the 12 punch will be the extra punch, and it 
does not correspond to any coefficient. 



SEMINAR PROCEEDINGS 
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Example 

The following determinant arose in the evaluation of 
the integral: 

7 

-14 -2i 



^ = 2- 



fc/ — 00 



D 



where 



I ''^■' 



dx 



,i=0 / Vi=0 



Ao A2 A^ Aq 

Ao A2 Ai Ae, 

Ao A2 ^4 Aq 

Ai As A5 Aj 

Ai A3 A5 Aj 

A^ As A5 A-! 



Ar = 



r 

2^ a^x^-"-' 



If the above evaluation method were applied to the determi- 
nant in this form, approximately 130,000 term cards and 



flo 


«2 + 0,1^ 


tti 


+ 


azx 


fle + cL^x 











«0 


a2 


+ 


axx 


a4 + asx 


Ao 











ao 






a2 + aix 


A, 


A 


«i 


a^ — ttxx^ 


fl5 


— 


asx" 


a? — a^x^ 











Ol 


as 


— 


a^x^ 


as — asx^ 


a-t 











ax 






as — aix^ 


as 


a^ 



24,576 selection cards would be required. However, by 
judicious manipulation of rows and columns, D can be 
written in the form: 



D = 



This form of D was evaluated by the above procedure. 
Notice that the zero elements are the same in both forms 
but that the degree of the elements appearing in the first 
four columns has been reduced. In this form, approximately 
3,000 term cards and 576 selection cards were required. 
The negative signs in columns 2, 3, and 4 were handled 
quite easily in the master gang punching step by double 
punching an X in the corresponding position of the term 
number. An extra crossfooting operation was used to de- 
termine the correct sign of the term card. The final answer 
involved about 340 terms associated with 19 different pow- 
ers of X. The entire problem required about 13 hours of 
machine time. 



